stanEx2
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}regression
ex5.stan
multi normal regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~normal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=normal_rng(m1[i],s);
}
}normal linear models
n=30
mdl=cmdstan_model('./ex5.stan')single regression
tb=tibble(x=runif(n,0,9),
y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)
qplot(data=tb,x,y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -798.475
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -14.6353 6.65847e-05 0.000885314 0.8846 0.8846 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -14.64
## b[1] 0.83
## b[2] 0.82
## s 0.99
## y1[1] 5.89
## y1[2] 4.33
## y1[3] 2.59
## y1[4] 2.39
## y1[5] 3.49
## y1[6] 5.18
## y1[7] 6.10
## y1[8] 2.58
## y1[9] 7.74
## y1[10] 3.37
## y1[11] 2.70
## y1[12] 2.57
## y1[13] 3.50
## y1[14] 0.49
## y1[15] 7.64
## y1[16] 7.02
## y1[17] 0.22
## y1[18] 6.68
## y1[19] 4.37
## y1[20] 5.99
## y1[21] 1.46
## y1[22] 7.57
## y1[23] 5.39
## y1[24] 6.89
## y1[25] 6.96
## y1[26] 4.34
## y1[27] 3.90
## y1[28] 8.82
## y1[29] 3.70
## y1[30] 6.12
## m1[1] 5.57
## m1[2] 4.32
## m1[3] 1.75
## m1[4] 2.84
## m1[5] 4.06
## m1[6] 4.33
## m1[7] 5.36
## m1[8] 1.68
## m1[9] 6.99
## m1[10] 3.06
## m1[11] 2.84
## m1[12] 3.55
## m1[13] 4.85
## m1[14] 1.74
## m1[15] 7.21
## m1[16] 7.95
## m1[17] 1.72
## m1[18] 8.00
## m1[19] 2.53
## m1[20] 6.36
## m1[21] 3.30
## m1[22] 7.04
## m1[23] 6.98
## m1[24] 6.44
## m1[25] 6.08
## m1[26] 4.66
## m1[27] 5.35
## m1[28] 8.02
## m1[29] 3.37
## m1[30] 6.67mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.17 -15.84 1.23 1.04 -18.63 -14.84 1.01 539 1012
## b[1] 0.83 0.83 0.42 0.41 0.13 1.51 1.01 520 527
## b[2] 0.82 0.82 0.08 0.08 0.69 0.95 1.00 623 773
## s 1.07 1.05 0.15 0.14 0.86 1.34 1.00 1290 1176
## y1[1] 5.56 5.55 1.07 1.05 3.84 7.30 1.00 2074 1969
## y1[2] 4.31 4.30 1.08 1.04 2.54 6.06 1.00 1853 1931
## y1[3] 1.75 1.75 1.13 1.09 -0.19 3.60 1.00 1814 1781
## y1[4] 2.80 2.79 1.11 1.09 0.95 4.63 1.00 2038 1768
## y1[5] 4.07 4.06 1.11 1.08 2.25 5.86 1.00 1808 1894
## y1[6] 4.32 4.33 1.11 1.03 2.45 6.13 1.00 1831 1828
## y1[7] 5.35 5.35 1.10 1.08 3.54 7.16 1.00 1924 1996
## y1[8] 1.72 1.72 1.12 1.12 -0.16 3.52 1.00 1463 1757
## y1[9] 7.00 6.98 1.09 1.07 5.23 8.77 1.00 2117 2015
## y1[10] 3.08 3.08 1.13 1.05 1.22 4.93 1.00 1945 1866
## y1[11] 2.87 2.86 1.11 1.07 1.07 4.74 1.00 1947 1846
## y1[12] 3.56 3.54 1.11 1.07 1.73 5.40 1.00 1900 1945
## y1[13] 4.87 4.87 1.10 1.11 3.07 6.63 1.00 2030 1731
## y1[14] 1.74 1.71 1.13 1.13 -0.08 3.58 1.00 1558 1769
## y1[15] 7.21 7.20 1.10 1.12 5.41 9.04 1.00 1958 1929
## y1[16] 7.96 7.96 1.16 1.13 6.06 9.84 1.00 1607 1601
## y1[17] 1.72 1.71 1.14 1.11 -0.10 3.59 1.00 1646 1840
## y1[18] 7.99 7.99 1.13 1.09 6.13 9.89 1.00 2186 1881
## y1[19] 2.50 2.49 1.10 1.07 0.66 4.34 1.00 1537 1897
## y1[20] 6.36 6.38 1.10 1.07 4.50 8.13 1.00 2212 1859
## y1[21] 3.27 3.30 1.12 1.06 1.43 5.08 1.00 1955 1721
## y1[22] 7.03 7.03 1.14 1.17 5.22 8.93 1.00 1940 1772
## y1[23] 6.99 6.98 1.13 1.10 5.18 8.78 1.00 2043 1858
## y1[24] 6.42 6.40 1.13 1.09 4.59 8.30 1.00 1751 1971
## y1[25] 6.14 6.13 1.12 1.10 4.33 7.97 1.00 2187 1930
## y1[26] 4.66 4.65 1.10 1.10 2.85 6.50 1.00 1917 1770
## y1[27] 5.36 5.39 1.11 1.10 3.53 7.10 1.00 1900 1772
## y1[28] 8.02 8.00 1.15 1.11 6.14 9.89 1.00 1962 1951
## y1[29] 3.40 3.41 1.09 1.05 1.67 5.19 1.00 1900 1828
## y1[30] 6.70 6.70 1.08 1.04 4.92 8.42 1.00 1987 1815
## m1[1] 5.57 5.57 0.21 0.21 5.23 5.92 1.00 2126 1452
## m1[2] 4.32 4.32 0.20 0.20 3.99 4.65 1.00 953 1037
## m1[3] 1.76 1.76 0.35 0.34 1.19 2.31 1.01 527 534
## m1[4] 2.84 2.84 0.27 0.27 2.40 3.27 1.01 570 705
## m1[5] 4.07 4.07 0.21 0.20 3.73 4.41 1.01 824 941
## m1[6] 4.33 4.34 0.20 0.20 4.01 4.67 1.00 963 1123
## m1[7] 5.37 5.37 0.20 0.20 5.03 5.70 1.00 2033 1474
## m1[8] 1.68 1.69 0.35 0.35 1.11 2.25 1.01 526 534
## m1[9] 6.99 6.99 0.29 0.28 6.52 7.47 1.00 1758 1408
## m1[10] 3.06 3.06 0.25 0.26 2.64 3.47 1.01 590 770
## m1[11] 2.84 2.84 0.27 0.27 2.40 3.28 1.01 570 705
## m1[12] 3.55 3.56 0.23 0.23 3.17 3.93 1.01 665 915
## m1[13] 4.85 4.85 0.20 0.20 4.53 5.18 1.00 1428 1397
## m1[14] 1.74 1.74 0.35 0.34 1.17 2.30 1.01 527 531
## m1[15] 7.21 7.21 0.30 0.30 6.73 7.72 1.00 1630 1398
## m1[16] 7.96 7.95 0.36 0.36 7.37 8.55 1.00 1300 1365
## m1[17] 1.73 1.73 0.35 0.34 1.16 2.29 1.01 527 534
## m1[18] 8.00 8.00 0.36 0.36 7.41 8.60 1.00 1287 1365
## m1[19] 2.53 2.53 0.29 0.29 2.05 2.99 1.01 551 641
## m1[20] 6.36 6.36 0.25 0.24 5.96 6.77 1.01 2038 1428
## m1[21] 3.30 3.30 0.24 0.24 2.90 3.69 1.01 621 893
## m1[22] 7.04 7.04 0.29 0.29 6.57 7.53 1.00 1730 1388
## m1[23] 6.98 6.98 0.29 0.28 6.51 7.46 1.00 1765 1408
## m1[24] 6.45 6.45 0.25 0.24 6.04 6.87 1.01 2010 1427
## m1[25] 6.08 6.08 0.23 0.23 5.70 6.47 1.01 2117 1400
## m1[26] 4.66 4.66 0.20 0.20 4.34 4.99 1.00 1219 1313
## m1[27] 5.35 5.35 0.20 0.20 5.01 5.68 1.00 2022 1473
## m1[28] 8.02 8.02 0.36 0.36 7.43 8.62 1.00 1282 1365
## m1[29] 3.37 3.37 0.24 0.24 2.97 3.76 1.01 632 901
## m1[30] 6.67 6.67 0.27 0.26 6.24 7.12 1.00 1920 1348
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
multiple regression
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -111.072
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 31 -9.1604 1.6178e-05 0.000692474 1 1 42
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.16
## b[1] -0.18
## b[2] 1.06
## b[3] -0.99
## s 0.82
## y1[1] 2.99
## y1[2] -1.19
## y1[3] 4.02
## y1[4] -5.16
## y1[5] 5.39
## y1[6] -1.92
## y1[7] -0.87
## y1[8] -2.35
## y1[9] -1.63
## y1[10] -3.67
## y1[11] -7.58
## y1[12] -0.71
## y1[13] -0.74
## y1[14] -0.21
## y1[15] -2.75
## y1[16] -6.00
## y1[17] -0.03
## y1[18] 0.40
## y1[19] 4.35
## y1[20] 0.18
## y1[21] -6.89
## y1[22] -2.82
## y1[23] 2.08
## y1[24] 0.14
## y1[25] -5.08
## y1[26] -4.82
## y1[27] -2.77
## y1[28] 0.38
## y1[29] 5.09
## y1[30] 4.05
## m1[1] 3.22
## m1[2] -0.57
## m1[3] 4.13
## m1[4] -4.61
## m1[5] 5.12
## m1[6] -0.99
## m1[7] -0.08
## m1[8] -2.43
## m1[9] -1.14
## m1[10] -3.57
## m1[11] -6.95
## m1[12] -1.69
## m1[13] -0.78
## m1[14] 0.33
## m1[15] -3.46
## m1[16] -7.04
## m1[17] 0.05
## m1[18] -0.70
## m1[19] 3.77
## m1[20] -0.60
## m1[21] -6.09
## m1[22] -3.89
## m1[23] 1.19
## m1[24] 0.43
## m1[25] -5.64
## m1[26] -4.16
## m1[27] -0.64
## m1[28] -1.17
## m1[29] 5.12
## m1[30] 3.67mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.38 -11.08 1.39 1.22 -14.05 -9.72 1.00 738 1115
## b[1] -0.20 -0.20 0.42 0.41 -0.89 0.48 1.01 382 413
## b[2] 1.06 1.06 0.06 0.06 0.95 1.16 1.00 815 1217
## b[3] -0.99 -0.99 0.06 0.06 -1.09 -0.88 1.00 682 1101
## s 0.91 0.90 0.12 0.12 0.73 1.14 1.00 1021 964
## y1[1] 3.25 3.24 0.96 0.93 1.68 4.82 1.00 2206 2039
## y1[2] -0.56 -0.57 0.95 0.95 -2.11 1.01 1.00 1706 1845
## y1[3] 4.12 4.13 0.96 0.93 2.54 5.66 1.00 1922 1745
## y1[4] -4.59 -4.60 0.96 0.94 -6.13 -3.01 1.00 1677 1599
## y1[5] 5.12 5.13 0.98 0.96 3.52 6.70 1.00 1760 1783
## y1[6] -1.03 -1.03 1.01 0.99 -2.67 0.65 1.00 1260 1688
## y1[7] -0.08 -0.07 0.98 0.99 -1.71 1.47 1.00 1705 1820
## y1[8] -2.40 -2.43 0.92 0.83 -3.86 -0.86 1.00 1880 1746
## y1[9] -1.13 -1.11 0.95 0.88 -2.71 0.43 1.00 1963 1842
## y1[10] -3.57 -3.57 0.96 0.92 -5.19 -2.04 1.00 1967 1994
## y1[11] -6.95 -6.97 0.97 0.93 -8.50 -5.33 1.00 2087 1769
## y1[12] -1.69 -1.66 1.00 0.95 -3.42 -0.10 1.00 1719 2016
## y1[13] -0.80 -0.81 0.96 0.95 -2.43 0.75 1.00 2110 1831
## y1[14] 0.30 0.28 0.97 0.95 -1.23 1.94 1.00 1858 1973
## y1[15] -3.46 -3.49 0.94 0.91 -4.96 -1.88 1.00 1890 1893
## y1[16] -7.04 -7.03 0.97 0.96 -8.61 -5.42 1.00 1693 1810
## y1[17] 0.07 0.05 0.93 0.93 -1.45 1.51 1.00 1935 1673
## y1[18] -0.74 -0.73 0.93 0.92 -2.20 0.78 1.00 1961 1774
## y1[19] 3.76 3.80 0.97 0.95 2.11 5.31 1.00 1909 1819
## y1[20] -0.60 -0.58 0.95 0.92 -2.15 0.97 1.00 1812 1892
## y1[21] -6.08 -6.07 0.98 0.99 -7.70 -4.53 1.00 1946 1972
## y1[22] -3.90 -3.89 0.97 0.97 -5.54 -2.33 1.00 1982 1855
## y1[23] 1.17 1.14 0.94 0.92 -0.40 2.69 1.00 1513 1874
## y1[24] 0.38 0.38 1.01 1.00 -1.28 2.11 1.00 1299 1711
## y1[25] -5.65 -5.63 0.96 0.93 -7.20 -4.09 1.00 1911 1951
## y1[26] -4.14 -4.11 0.96 0.93 -5.74 -2.58 1.00 1975 1659
## y1[27] -0.62 -0.61 0.97 0.96 -2.26 0.89 1.00 1684 2033
## y1[28] -1.19 -1.19 0.92 0.90 -2.69 0.28 1.00 1915 1854
## y1[29] 5.09 5.09 1.01 0.96 3.47 6.78 1.00 1952 2042
## y1[30] 3.69 3.70 0.95 0.93 2.11 5.22 1.00 1724 1802
## m1[1] 3.23 3.24 0.29 0.28 2.75 3.71 1.00 1819 1619
## m1[2] -0.58 -0.58 0.28 0.28 -1.05 -0.11 1.01 401 520
## m1[3] 4.14 4.15 0.33 0.32 3.59 4.68 1.00 1760 1563
## m1[4] -4.62 -4.62 0.27 0.26 -5.07 -4.18 1.01 700 1159
## m1[5] 5.12 5.11 0.35 0.34 4.56 5.69 1.00 1199 1401
## m1[6] -1.01 -1.01 0.35 0.34 -1.58 -0.43 1.01 385 445
## m1[7] -0.09 -0.09 0.27 0.27 -0.52 0.35 1.01 418 547
## m1[8] -2.44 -2.44 0.19 0.20 -2.77 -2.14 1.01 722 1069
## m1[9] -1.13 -1.13 0.22 0.22 -1.50 -0.76 1.00 1705 1385
## m1[10] -3.57 -3.57 0.24 0.23 -3.95 -3.16 1.00 1941 1740
## m1[11] -6.95 -6.94 0.34 0.33 -7.50 -6.40 1.00 1822 1695
## m1[12] -1.68 -1.68 0.31 0.30 -2.18 -1.17 1.00 870 1286
## m1[13] -0.78 -0.77 0.24 0.24 -1.18 -0.37 1.00 1280 1347
## m1[14] 0.33 0.33 0.22 0.21 -0.03 0.70 1.00 1840 1596
## m1[15] -3.45 -3.46 0.25 0.24 -3.85 -3.03 1.00 1651 1789
## m1[16] -7.04 -7.03 0.33 0.33 -7.59 -6.50 1.00 1742 1511
## m1[17] 0.06 0.07 0.26 0.26 -0.38 0.51 1.00 1151 1595
## m1[18] -0.70 -0.69 0.17 0.16 -0.97 -0.42 1.00 1693 1420
## m1[19] 3.78 3.79 0.34 0.33 3.21 4.33 1.00 1426 1545
## m1[20] -0.61 -0.61 0.24 0.25 -1.01 -0.21 1.01 432 662
## m1[21] -6.10 -6.09 0.31 0.29 -6.61 -5.61 1.00 1042 1255
## m1[22] -3.89 -3.89 0.22 0.22 -4.24 -3.53 1.00 1968 1639
## m1[23] 1.18 1.19 0.28 0.28 0.73 1.63 1.01 461 718
## m1[24] 0.42 0.42 0.38 0.37 -0.21 1.04 1.01 389 394
## m1[25] -5.64 -5.64 0.28 0.27 -6.10 -5.19 1.00 1748 1438
## m1[26] -4.16 -4.17 0.22 0.22 -4.53 -3.80 1.00 1437 1484
## m1[27] -0.63 -0.62 0.34 0.33 -1.21 -0.07 1.00 772 1253
## m1[28] -1.17 -1.17 0.17 0.17 -1.46 -0.90 1.00 953 1197
## m1[29] 5.12 5.12 0.35 0.35 4.56 5.70 1.00 1664 1420
## m1[30] 3.67 3.66 0.29 0.29 3.18 4.15 1.00 926 1380
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANOVA
tb=tibble(c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)
qplot(data=tb,c,y,geom='boxplot')
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -311.803
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -12.4643 5.40671e-05 0.000515921 0.9009 0.9009 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.46
## b[1] 0.15
## b[2] 2.05
## b[3] -2.02
## s 0.92
## y1[1] 3.55
## y1[2] -1.81
## y1[3] 0.17
## y1[4] 0.73
## y1[5] 1.75
## y1[6] 1.00
## y1[7] -1.14
## y1[8] 0.64
## y1[9] -0.77
## y1[10] 1.05
## y1[11] 0.94
## y1[12] 1.39
## y1[13] 0.22
## y1[14] 0.62
## y1[15] 1.26
## y1[16] 1.60
## y1[17] 1.74
## y1[18] 2.78
## y1[19] -0.82
## y1[20] 1.84
## y1[21] -2.87
## y1[22] 1.70
## y1[23] -2.17
## y1[24] -2.73
## y1[25] 1.83
## y1[26] -1.98
## y1[27] 2.11
## y1[28] 1.48
## y1[29] 2.02
## y1[30] -0.82
## m1[1] 2.19
## m1[2] -1.87
## m1[3] 0.15
## m1[4] 2.19
## m1[5] 0.15
## m1[6] 2.19
## m1[7] 0.15
## m1[8] 0.15
## m1[9] 0.15
## m1[10] 0.15
## m1[11] 0.15
## m1[12] 0.15
## m1[13] 0.15
## m1[14] 0.15
## m1[15] 0.15
## m1[16] 2.19
## m1[17] 2.19
## m1[18] 2.19
## m1[19] 0.15
## m1[20] 2.19
## m1[21] -1.87
## m1[22] 2.19
## m1[23] -1.87
## m1[24] -1.87
## m1[25] 2.19
## m1[26] -1.87
## m1[27] 2.19
## m1[28] 2.19
## m1[29] 2.19
## m1[30] -1.87mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.78 -14.46 1.52 1.33 -17.79 -12.95 1.00 736 992
## b[1] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## b[2] 2.04 2.02 0.42 0.43 1.36 2.75 1.00 927 910
## b[3] -2.02 -2.01 0.52 0.52 -2.88 -1.18 1.00 762 945
## s 1.02 1.00 0.15 0.14 0.80 1.29 1.00 2269 1757
## y1[1] 2.18 2.20 1.09 1.07 0.37 3.95 1.00 1988 1906
## y1[2] -1.84 -1.85 1.11 1.07 -3.73 -0.03 1.00 1606 1894
## y1[3] 0.17 0.19 1.09 1.08 -1.57 1.94 1.00 2122 2152
## y1[4] 2.18 2.20 1.09 1.07 0.38 3.94 1.00 1944 1959
## y1[5] 0.15 0.12 1.05 1.01 -1.57 1.90 1.00 1741 1980
## y1[6] 2.23 2.22 1.07 1.03 0.52 3.99 1.00 1796 1852
## y1[7] 0.17 0.16 1.07 1.02 -1.64 1.93 1.00 2066 1894
## y1[8] 0.16 0.13 1.08 1.02 -1.58 1.96 1.00 1931 1847
## y1[9] 0.13 0.11 1.06 1.04 -1.58 1.91 1.00 1719 1803
## y1[10] 0.17 0.18 1.06 1.05 -1.60 1.84 1.00 1883 1715
## y1[11] 0.13 0.11 1.08 1.08 -1.62 1.90 1.00 1943 1877
## y1[12] 0.18 0.17 1.08 1.05 -1.53 2.03 1.00 1989 1811
## y1[13] 0.19 0.18 1.08 1.09 -1.56 1.90 1.00 1890 1888
## y1[14] 0.19 0.21 1.05 1.02 -1.59 1.86 1.00 1808 1934
## y1[15] 0.15 0.15 1.07 1.02 -1.65 1.90 1.00 1776 1933
## y1[16] 2.19 2.18 1.07 1.07 0.44 3.96 1.00 1841 1945
## y1[17] 2.18 2.20 1.10 1.09 0.37 3.90 1.00 1973 1732
## y1[18] 2.18 2.13 1.09 1.03 0.40 3.98 1.00 1698 1820
## y1[19] 0.13 0.15 1.07 1.03 -1.68 1.86 1.00 1950 1885
## y1[20] 2.23 2.23 1.05 1.02 0.51 3.92 1.00 1962 1781
## y1[21] -1.89 -1.91 1.14 1.08 -3.78 -0.01 1.00 1840 1915
## y1[22] 2.18 2.21 1.09 1.07 0.34 3.87 1.00 1932 1934
## y1[23] -1.89 -1.90 1.09 1.06 -3.70 -0.08 1.00 1715 1942
## y1[24] -1.83 -1.82 1.09 1.05 -3.65 -0.01 1.00 1918 1981
## y1[25] 2.22 2.22 1.08 1.07 0.53 4.00 1.00 1725 1878
## y1[26] -1.85 -1.84 1.12 1.09 -3.72 0.01 1.00 1882 1676
## y1[27] 2.20 2.20 1.03 1.02 0.55 3.93 1.00 1990 1894
## y1[28] 2.23 2.23 1.06 1.05 0.49 3.94 1.00 1974 1853
## y1[29] 2.14 2.14 1.10 1.06 0.36 4.02 1.00 2038 2040
## y1[30] -1.89 -1.89 1.14 1.11 -3.68 -0.04 1.00 1832 1786
## m1[1] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[2] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
## m1[3] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[4] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[5] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[6] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[7] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[8] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[9] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[10] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[11] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[12] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[13] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[14] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[15] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[16] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[17] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[18] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[19] 0.16 0.17 0.29 0.27 -0.32 0.65 1.00 980 1057
## m1[20] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[21] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
## m1[22] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[23] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
## m1[24] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
## m1[25] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[26] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
## m1[27] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[28] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[29] 2.20 2.19 0.31 0.31 1.70 2.71 1.00 1410 1404
## m1[30] -1.86 -1.85 0.44 0.45 -2.57 -1.13 1.00 1004 924
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANCOVA
tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))
f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)
lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])
qplot(data=tb,x,y,shape=c,size=I(2))
plot(tb$x,tb$y,col=1+lv[factor(tb$c)])
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -44312
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 42 -14.4242 0.000147739 0.00233842 0.6781 0.6781 94
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.42
## b[1] -0.04
## b[2] 1.06
## b[3] 1.38
## b[4] -2.43
## s 0.98
## y1[1] 5.22
## y1[2] 0.54
## y1[3] 0.88
## y1[4] 4.00
## y1[5] 4.53
## y1[6] 1.56
## y1[7] 9.10
## y1[8] 4.33
## y1[9] -0.57
## y1[10] 1.62
## y1[11] 5.38
## y1[12] 0.29
## y1[13] 5.83
## y1[14] 2.24
## y1[15] 5.13
## y1[16] -0.73
## y1[17] 2.41
## y1[18] 8.63
## y1[19] 6.02
## y1[20] 4.08
## y1[21] 3.69
## y1[22] 7.15
## y1[23] -4.73
## y1[24] 6.66
## y1[25] 8.62
## y1[26] 6.05
## y1[27] 10.11
## y1[28] 5.76
## y1[29] 1.72
## y1[30] 5.00
## m1[1] 5.69
## m1[2] 1.68
## m1[3] 0.57
## m1[4] 3.37
## m1[5] 3.89
## m1[6] -0.14
## m1[7] 9.54
## m1[8] 4.88
## m1[9] -0.37
## m1[10] 0.21
## m1[11] 6.47
## m1[12] 0.87
## m1[13] 6.61
## m1[14] 2.10
## m1[15] 5.25
## m1[16] 0.42
## m1[17] 1.24
## m1[18] 9.10
## m1[19] 4.96
## m1[20] 4.63
## m1[21] 2.30
## m1[22] 6.84
## m1[23] -2.07
## m1[24] 7.42
## m1[25] 7.68
## m1[26] 7.63
## m1[27] 8.21
## m1[28] 4.64
## m1[29] 2.54
## m1[30] 4.95mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.35 -16.93 1.89 1.61 -20.93 -15.08 1.00 713 855
## b[1] -0.03 -0.03 0.50 0.47 -0.87 0.83 1.00 680 805
## b[2] 1.05 1.05 0.09 0.09 0.91 1.19 1.00 1166 1111
## b[3] 1.39 1.40 0.53 0.53 0.52 2.29 1.01 667 350
## b[4] -2.43 -2.45 0.48 0.46 -3.21 -1.65 1.00 714 919
## s 1.12 1.10 0.17 0.16 0.88 1.42 1.00 1496 1120
## y1[1] 5.68 5.69 1.22 1.23 3.68 7.64 1.00 1857 1864
## y1[2] 1.73 1.72 1.17 1.15 -0.19 3.60 1.00 1689 1742
## y1[3] 0.62 0.63 1.20 1.17 -1.34 2.54 1.00 1579 1847
## y1[4] 3.42 3.42 1.24 1.20 1.34 5.47 1.00 1798 1768
## y1[5] 3.84 3.85 1.17 1.15 1.89 5.72 1.00 1853 1856
## y1[6] -0.14 -0.15 1.18 1.15 -2.14 1.70 1.00 1884 1761
## y1[7] 9.55 9.53 1.27 1.27 7.48 11.61 1.00 1819 1784
## y1[8] 4.85 4.87 1.22 1.15 2.86 6.84 1.00 1809 1819
## y1[9] -0.38 -0.41 1.20 1.13 -2.33 1.61 1.00 1817 2050
## y1[10] 0.20 0.19 1.14 1.13 -1.66 1.99 1.00 1612 1928
## y1[11] 6.45 6.45 1.21 1.16 4.44 8.40 1.00 1957 2001
## y1[12] 0.84 0.80 1.21 1.13 -1.13 2.89 1.00 2030 1775
## y1[13] 6.62 6.64 1.18 1.08 4.64 8.52 1.00 1906 1931
## y1[14] 2.07 2.07 1.19 1.16 0.14 4.04 1.00 1972 1904
## y1[15] 5.20 5.21 1.20 1.12 3.18 7.12 1.00 1928 1859
## y1[16] 0.41 0.40 1.21 1.19 -1.56 2.41 1.00 1826 1944
## y1[17] 1.23 1.25 1.21 1.20 -0.76 3.21 1.00 1647 1933
## y1[18] 9.11 9.13 1.21 1.19 7.13 11.08 1.00 1945 1839
## y1[19] 4.98 5.00 1.17 1.10 3.01 6.86 1.00 2071 1807
## y1[20] 4.61 4.59 1.18 1.15 2.73 6.57 1.00 1857 1708
## y1[21] 2.31 2.31 1.19 1.18 0.39 4.20 1.00 1715 1704
## y1[22] 6.84 6.85 1.20 1.15 4.99 8.83 1.00 1862 1967
## y1[23] -2.05 -2.07 1.23 1.20 -4.03 0.01 1.00 1855 1674
## y1[24] 7.42 7.44 1.24 1.21 5.43 9.38 1.00 1814 1806
## y1[25] 7.68 7.67 1.21 1.13 5.71 9.72 1.00 1602 1944
## y1[26] 7.63 7.61 1.21 1.18 5.60 9.56 1.00 1936 1816
## y1[27] 8.26 8.26 1.24 1.20 6.15 10.27 1.00 1328 1258
## y1[28] 4.64 4.67 1.18 1.16 2.72 6.53 1.00 1768 1848
## y1[29] 2.50 2.50 1.26 1.20 0.49 4.56 1.00 1743 1533
## y1[30] 4.93 4.95 1.25 1.19 2.88 6.93 1.00 2054 1981
## m1[1] 5.68 5.69 0.48 0.45 4.91 6.46 1.00 1180 1539
## m1[2] 1.69 1.69 0.41 0.38 1.00 2.37 1.00 644 770
## m1[3] 0.58 0.59 0.47 0.43 -0.21 1.37 1.00 664 829
## m1[4] 3.39 3.39 0.50 0.49 2.56 4.21 1.00 1029 832
## m1[5] 3.89 3.90 0.34 0.32 3.31 4.43 1.00 688 713
## m1[6] -0.14 -0.13 0.39 0.36 -0.79 0.51 1.00 1236 1344
## m1[7] 9.54 9.53 0.51 0.48 8.70 10.39 1.00 1098 1159
## m1[8] 4.86 4.88 0.44 0.42 4.15 5.58 1.00 1193 1552
## m1[9] -0.37 -0.37 0.40 0.37 -1.04 0.28 1.00 1242 1343
## m1[10] 0.22 0.22 0.38 0.35 -0.42 0.84 1.00 1228 1308
## m1[11] 6.48 6.48 0.44 0.42 5.73 7.19 1.00 1036 821
## m1[12] 0.87 0.87 0.36 0.34 0.26 1.46 1.00 1221 1326
## m1[13] 6.60 6.61 0.38 0.36 5.98 7.22 1.00 1023 1096
## m1[14] 2.09 2.10 0.36 0.34 1.50 2.67 1.00 1245 1257
## m1[15] 5.23 5.25 0.46 0.43 4.50 5.98 1.00 1191 1531
## m1[16] 0.42 0.42 0.37 0.34 -0.20 1.03 1.00 1218 1364
## m1[17] 1.25 1.25 0.43 0.40 0.52 1.97 1.00 650 823
## m1[18] 9.08 9.08 0.50 0.48 8.28 9.92 1.00 1342 1137
## m1[19] 4.95 4.96 0.34 0.32 4.36 5.49 1.00 783 856
## m1[20] 4.63 4.63 0.34 0.32 4.05 5.17 1.00 748 782
## m1[21] 2.31 2.31 0.38 0.36 1.67 2.94 1.00 641 711
## m1[22] 6.83 6.84 0.39 0.38 6.20 7.46 1.00 1061 1188
## m1[23] -2.06 -2.07 0.47 0.45 -2.83 -1.30 1.00 1213 1175
## m1[24] 7.41 7.41 0.41 0.41 6.75 8.08 1.00 1148 1318
## m1[25] 7.69 7.68 0.45 0.43 6.93 8.43 1.00 1056 848
## m1[26] 7.63 7.63 0.45 0.43 6.88 8.37 1.00 1055 839
## m1[27] 8.22 8.21 0.46 0.44 7.44 8.98 1.00 1070 864
## m1[28] 4.63 4.64 0.34 0.32 4.05 5.17 1.00 749 782
## m1[29] 2.56 2.57 0.53 0.52 1.68 3.46 1.00 1039 940
## m1[30] 4.93 4.95 0.44 0.42 4.22 5.65 1.00 1192 1532
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
generalized linear regression
poisson regression
objective variable [0,Infinity), integer. variance of error is near to mean
(normal linear regression, correlation is near to 1,-1, variance of error is constant)
# of samples is N,
log li=sum(bj*xji),j=[0,K],i=[1,N]
yi~Po(li)
or
li=sum(bj xji),j=[0,k]
yi~Po(exp li)
when xj -> xj +1, y -> y* exp bj ex6-1.stan
poisson regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] l=X*b;
y~poisson_log(l);
}
generated quantities{
array[N] int y1;
vector[N] l1=X*b;
for(i in 1:N){
y1[i]=poisson_log_rng(l1[i]);
}
}n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='poisson')##
## Call: glm(formula = f, family = "poisson", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.035 0.822 0.344
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 36.7
## Residual Deviance: 27.3 AIC: 94.7X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -77.8209
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -19.272 0.000160908 0.000572024 0.8673 0.8673 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -19.27
## b[1] 0.03
## b[2] 0.82
## b[3] 0.34
## y1[1] 2.00
## y1[2] 1.00
## y1[3] 2.00
## y1[4] 2.00
## y1[5] 5.00
## y1[6] 0.00
## y1[7] 0.00
## y1[8] 0.00
## y1[9] 3.00
## y1[10] 1.00
## y1[11] 0.00
## y1[12] 1.00
## y1[13] 1.00
## y1[14] 2.00
## y1[15] 1.00
## y1[16] 0.00
## y1[17] 2.00
## y1[18] 0.00
## y1[19] 2.00
## y1[20] 0.00
## y1[21] 2.00
## y1[22] 1.00
## y1[23] 3.00
## y1[24] 0.00
## y1[25] 1.00
## y1[26] 5.00
## y1[27] 2.00
## y1[28] 2.00
## y1[29] 3.00
## y1[30] 4.00
## l1[1] 0.04
## l1[2] 0.37
## l1[3] 0.86
## l1[4] 0.90
## l1[5] 0.71
## l1[6] 0.70
## l1[7] 0.57
## l1[8] 0.59
## l1[9] 0.78
## l1[10] -0.65
## l1[11] -0.59
## l1[12] -0.77
## l1[13] 0.59
## l1[14] 0.50
## l1[15] 0.46
## l1[16] 0.58
## l1[17] -0.08
## l1[18] 0.15
## l1[19] 0.82
## l1[20] 0.52
## l1[21] -0.36
## l1[22] 0.50
## l1[23] 1.07
## l1[24] -0.20
## l1[25] 0.26
## l1[26] 0.65
## l1[27] 0.80
## l1[28] 0.70
## l1[29] 1.12
## l1[30] 1.08mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='l1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -20.81 -20.49 1.29 1.03 -23.42 -19.45 1.01 816 1016
## b[1] -0.01 0.01 0.29 0.28 -0.54 0.43 1.02 669 826
## b[2] 0.85 0.85 0.32 0.33 0.34 1.40 1.01 785 977
## b[3] 0.36 0.34 0.30 0.29 -0.11 0.86 1.01 796 912
## y1[1] 1.01 1.00 1.06 1.48 0.00 3.00 1.00 1892 1923
## y1[2] 1.42 1.00 1.22 1.48 0.00 4.00 1.00 2005 1958
## y1[3] 2.36 2.00 1.69 1.48 0.00 6.00 1.00 1796 1678
## y1[4] 2.47 2.00 1.70 1.48 0.00 6.00 1.00 2219 1913
## y1[5] 2.03 2.00 1.51 1.48 0.00 5.00 1.00 2057 1951
## y1[6] 2.03 2.00 1.50 1.48 0.00 5.00 1.00 1772 1824
## y1[7] 1.75 2.00 1.33 1.48 0.00 4.00 1.00 1961 1728
## y1[8] 1.84 2.00 1.39 1.48 0.00 4.00 1.00 1884 1770
## y1[9] 2.12 2.00 1.53 1.48 0.00 5.00 1.00 1868 2035
## y1[10] 0.54 0.00 0.78 0.00 0.00 2.00 1.00 1753 1847
## y1[11] 0.56 0.00 0.78 0.00 0.00 2.00 1.00 1917 1745
## y1[12] 0.47 0.00 0.72 0.00 0.00 2.00 1.00 1660 1694
## y1[13] 1.81 2.00 1.41 1.48 0.00 4.00 1.00 1680 1700
## y1[14] 1.66 1.00 1.34 1.48 0.00 4.00 1.00 1705 1801
## y1[15] 1.58 1.00 1.28 1.48 0.00 4.00 1.00 2003 1845
## y1[16] 1.72 2.00 1.32 1.48 0.00 4.00 1.00 1904 1904
## y1[17] 0.95 1.00 1.03 1.48 0.00 3.00 1.00 1860 1703
## y1[18] 1.17 1.00 1.17 1.48 0.00 3.00 1.00 1982 2021
## y1[19] 2.28 2.00 1.61 1.48 0.00 5.00 1.00 1749 1564
## y1[20] 1.70 1.00 1.37 1.48 0.00 4.00 1.00 1806 1915
## y1[21] 0.69 0.00 0.88 0.00 0.00 2.00 1.00 1852 1902
## y1[22] 1.61 1.00 1.34 1.48 0.00 4.00 1.00 1743 1931
## y1[23] 2.98 3.00 1.85 1.48 0.00 6.00 1.00 1670 1960
## y1[24] 0.85 1.00 0.99 1.48 0.00 3.00 1.00 1782 1936
## y1[25] 1.28 1.00 1.14 1.48 0.00 3.00 1.00 1842 1796
## y1[26] 1.96 2.00 1.48 1.48 0.00 5.00 1.00 2084 2058
## y1[27] 2.28 2.00 1.64 1.48 0.00 5.00 1.00 1770 2010
## y1[28] 2.00 2.00 1.50 1.48 0.00 5.00 1.00 1790 1856
## y1[29] 3.10 3.00 1.88 1.48 0.00 7.00 1.00 1674 1492
## y1[30] 2.96 3.00 1.84 1.48 0.00 6.00 1.00 1933 1985
## l1[1] -0.01 0.00 0.31 0.31 -0.54 0.48 1.01 1053 1124
## l1[2] 0.34 0.34 0.22 0.22 -0.03 0.69 1.00 1450 1474
## l1[3] 0.84 0.84 0.19 0.18 0.51 1.15 1.00 1468 1381
## l1[4] 0.89 0.89 0.19 0.18 0.55 1.21 1.00 1367 1407
## l1[5] 0.68 0.70 0.25 0.25 0.28 1.07 1.00 1203 1153
## l1[6] 0.68 0.68 0.18 0.17 0.39 0.97 1.00 1789 1387
## l1[7] 0.54 0.54 0.19 0.18 0.24 0.84 1.00 1781 1652
## l1[8] 0.57 0.58 0.24 0.23 0.18 0.93 1.01 1088 998
## l1[9] 0.76 0.76 0.18 0.17 0.45 1.06 1.00 1655 1408
## l1[10] -0.72 -0.69 0.50 0.49 -1.60 0.06 1.02 645 884
## l1[11] -0.66 -0.63 0.48 0.47 -1.51 0.08 1.02 642 877
## l1[12] -0.84 -0.82 0.54 0.53 -1.80 0.00 1.02 651 863
## l1[13] 0.57 0.58 0.24 0.23 0.18 0.93 1.01 1087 1009
## l1[14] 0.47 0.48 0.23 0.23 0.07 0.83 1.01 969 926
## l1[15] 0.42 0.43 0.20 0.20 0.08 0.75 1.00 1589 1434
## l1[16] 0.55 0.55 0.18 0.18 0.25 0.86 1.00 1792 1544
## l1[17] -0.13 -0.12 0.35 0.36 -0.72 0.42 1.01 998 983
## l1[18] 0.11 0.12 0.28 0.28 -0.37 0.55 1.01 1150 1147
## l1[19] 0.80 0.80 0.18 0.17 0.48 1.11 1.00 1558 1387
## l1[20] 0.49 0.50 0.23 0.23 0.09 0.85 1.01 994 947
## l1[21] -0.42 -0.40 0.41 0.40 -1.13 0.20 1.02 638 936
## l1[22] 0.47 0.49 0.23 0.23 0.08 0.84 1.01 976 945
## l1[23] 1.06 1.06 0.23 0.22 0.67 1.42 1.00 1099 1194
## l1[24] -0.25 -0.23 0.39 0.40 -0.91 0.36 1.01 960 1032
## l1[25] 0.22 0.23 0.25 0.25 -0.20 0.61 1.00 1285 1327
## l1[26] 0.63 0.65 0.24 0.24 0.23 1.00 1.00 1153 1024
## l1[27] 0.78 0.79 0.26 0.26 0.34 1.19 1.00 1235 1160
## l1[28] 0.68 0.69 0.25 0.25 0.27 1.06 1.00 1196 1153
## l1[29] 1.11 1.12 0.24 0.24 0.71 1.51 1.00 1042 1102
## l1[30] 1.07 1.07 0.23 0.23 0.68 1.44 1.00 1101 1193
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3
## 0 2 2 2 0
## 1 1 3 4 0
## 2 1 5 3 1
## 3 0 0 3 0
## 4 0 0 0 1
## 5 0 0 1 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3
## 0 2 0 1 0
## 1 1 0 2 0
## 2 1 3 1 0
## 3 0 0 2 0
## 4 0 0 0 0
## 5 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3
## 0 0 2 1 0
## 1 0 3 2 0
## 2 0 2 2 1
## 3 0 0 1 0
## 4 0 0 0 1
## 5 0 0 1 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3
## 0 3 3 0 0
## 1 4 3 1 0
## 2 1 7 1 1
## 3 0 1 2 0
## 4 0 0 0 1
## 5 0 0 2 0cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3
## 0 2 1 0 0
## 1 1 2 0 0
## 2 1 4 0 0
## 3 0 1 1 0
## 4 0 0 0 0
## 5 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3
## 0 1 2 0 0
## 1 3 1 1 0
## 2 0 3 1 1
## 3 0 0 1 0
## 4 0 0 0 1
## 5 0 0 2 0par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
logistic regression
# of samples is N,
objective variable 0/1 binary
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~Ber(pi), 0/1 binary
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-2.stan
logistic regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~bernoulli_logit(z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=bernoulli_rng(inv_logit(z1[i]));
}
}n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='binomial') # it can caluculte when all trials are once##
## Call: glm(formula = f, family = "binomial", data = tb)
##
## Coefficients:
## (Intercept) x cb
## -0.890 0.645 0.937
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 41.1
## Residual Deviance: 38.6 AIC: 44.6X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -37.8535
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -19.2851 0.000136794 3.65965e-05 0.9577 0.9577 15
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -19.29
## b[1] -0.89
## b[2] 0.64
## b[3] 0.94
## y1[1] 0.00
## y1[2] 0.00
## y1[3] 1.00
## y1[4] 0.00
## y1[5] 1.00
## y1[6] 0.00
## y1[7] 0.00
## y1[8] 1.00
## y1[9] 1.00
## y1[10] 1.00
## y1[11] 1.00
## y1[12] 1.00
## y1[13] 0.00
## y1[14] 0.00
## y1[15] 1.00
## y1[16] 0.00
## y1[17] 1.00
## y1[18] 0.00
## y1[19] 1.00
## y1[20] 0.00
## y1[21] 1.00
## y1[22] 0.00
## y1[23] 0.00
## y1[24] 0.00
## y1[25] 1.00
## y1[26] 1.00
## y1[27] 1.00
## y1[28] 0.00
## y1[29] 1.00
## y1[30] 0.00
## z1[1] -0.94
## z1[2] 0.54
## z1[3] -1.18
## z1[4] -0.83
## z1[5] 0.45
## z1[6] -0.33
## z1[7] -0.26
## z1[8] 0.49
## z1[9] -0.17
## z1[10] 0.13
## z1[11] -0.81
## z1[12] -0.31
## z1[13] 0.56
## z1[14] 0.35
## z1[15] 0.01
## z1[16] -1.43
## z1[17] 0.29
## z1[18] -0.03
## z1[19] 0.01
## z1[20] 0.67
## z1[21] 0.07
## z1[22] -0.57
## z1[23] 0.13
## z1[24] -1.20
## z1[25] -1.11
## z1[26] -0.59
## z1[27] -0.45
## z1[28] -0.56
## z1[29] -1.48
## z1[30] -0.37mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -20.86 -20.50 1.30 1.04 -23.43 -19.46 1.00 762 991
## b[1] -1.00 -0.96 0.67 0.66 -2.17 0.03 1.00 957 881
## b[2] 0.75 0.73 0.71 0.68 -0.38 1.96 1.00 1249 1041
## b[3] 1.03 1.00 0.82 0.80 -0.25 2.41 1.01 900 737
## y1[1] 0.26 0.00 0.44 0.00 0.00 1.00 1.00 1794 NA
## y1[2] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 1912 NA
## y1[3] 0.22 0.00 0.42 0.00 0.00 1.00 1.00 1718 NA
## y1[4] 0.30 0.00 0.46 0.00 0.00 1.00 1.00 1988 NA
## y1[5] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 1672 NA
## y1[6] 0.43 0.00 0.50 0.00 0.00 1.00 1.00 2053 NA
## y1[7] 0.44 0.00 0.50 0.00 0.00 1.00 1.00 1810 NA
## y1[8] 0.61 1.00 0.49 0.00 0.00 1.00 1.00 2153 NA
## y1[9] 0.43 0.00 0.50 0.00 0.00 1.00 1.00 1927 NA
## y1[10] 0.55 1.00 0.50 0.00 0.00 1.00 1.00 2113 NA
## y1[11] 0.31 0.00 0.46 0.00 0.00 1.00 1.00 2022 NA
## y1[12] 0.40 0.00 0.49 0.00 0.00 1.00 1.00 1811 NA
## y1[13] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 2035 NA
## y1[14] 0.60 1.00 0.49 0.00 0.00 1.00 1.00 1904 NA
## y1[15] 0.51 1.00 0.50 0.00 0.00 1.00 1.00 1949 NA
## y1[16] 0.21 0.00 0.41 0.00 0.00 1.00 1.00 1942 NA
## y1[17] 0.57 1.00 0.49 0.00 0.00 1.00 1.00 1908 NA
## y1[18] 0.48 0.00 0.50 0.00 0.00 1.00 1.00 1845 NA
## y1[19] 0.49 0.00 0.50 0.00 0.00 1.00 1.00 1938 NA
## y1[20] 0.64 1.00 0.48 0.00 0.00 1.00 1.00 2129 NA
## y1[21] 0.51 1.00 0.50 0.00 0.00 1.00 1.00 1741 NA
## y1[22] 0.34 0.00 0.47 0.00 0.00 1.00 1.00 1913 NA
## y1[23] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 2129 NA
## y1[24] 0.23 0.00 0.42 0.00 0.00 1.00 1.00 1893 NA
## y1[25] 0.25 0.00 0.43 0.00 0.00 1.00 1.00 1771 NA
## y1[26] 0.35 0.00 0.48 0.00 0.00 1.00 1.00 1735 NA
## y1[27] 0.40 0.00 0.49 0.00 0.00 1.00 1.00 1574 NA
## y1[28] 0.34 0.00 0.47 0.00 0.00 1.00 1.00 1938 NA
## y1[29] 0.20 0.00 0.40 0.00 0.00 1.00 1.00 1877 NA
## y1[30] 0.40 0.00 0.49 0.00 0.00 1.00 1.00 1903 NA
## z1[1] -1.07 -1.03 0.69 0.67 -2.30 -0.02 1.00 949 881
## z1[2] 0.59 0.56 0.69 0.66 -0.48 1.78 1.00 1591 1499
## z1[3] -1.35 -1.29 0.81 0.79 -2.77 -0.14 1.00 957 919
## z1[4] -0.93 -0.89 0.66 0.64 -2.06 0.08 1.00 971 934
## z1[5] 0.49 0.47 0.63 0.60 -0.49 1.58 1.00 1687 1404
## z1[6] -0.35 -0.34 0.79 0.79 -1.68 0.93 1.00 1233 1218
## z1[7] -0.27 -0.25 0.84 0.82 -1.69 1.11 1.00 1269 1219
## z1[8] 0.54 0.52 0.66 0.63 -0.47 1.68 1.00 1635 1414
## z1[9] -0.23 -0.22 0.61 0.59 -1.27 0.78 1.00 1830 1392
## z1[10] 0.12 0.11 0.52 0.51 -0.76 0.98 1.00 2093 1166
## z1[11] -0.91 -0.87 0.66 0.64 -2.05 0.10 1.00 975 916
## z1[12] -0.39 -0.38 0.70 0.68 -1.57 0.77 1.00 1677 1310
## z1[13] 0.63 0.59 0.71 0.68 -0.48 1.86 1.00 1567 1573
## z1[14] 0.38 0.36 0.58 0.54 -0.55 1.39 1.00 1823 1269
## z1[15] -0.02 -0.02 0.53 0.53 -0.92 0.87 1.00 2049 1127
## z1[16] -1.63 -1.58 0.99 0.96 -3.34 -0.16 1.00 990 891
## z1[17] 0.30 0.29 0.55 0.53 -0.59 1.25 1.00 1923 1176
## z1[18] -0.07 -0.06 0.55 0.54 -1.00 0.84 1.00 2013 1263
## z1[19] -0.02 -0.02 0.53 0.53 -0.92 0.87 1.00 2052 1127
## z1[20] 0.75 0.72 0.79 0.76 -0.50 2.11 1.00 1490 1683
## z1[21] 0.05 0.04 0.52 0.52 -0.86 0.91 1.00 2094 1266
## z1[22] -0.63 -0.60 0.68 0.66 -1.81 0.42 1.00 1093 1086
## z1[23] 0.12 0.12 0.52 0.51 -0.76 0.99 1.00 2091 1164
## z1[24] -1.36 -1.31 0.82 0.79 -2.82 -0.13 1.00 958 925
## z1[25] -1.26 -1.21 0.77 0.75 -2.61 -0.10 1.00 950 919
## z1[26] -0.72 -0.70 0.94 0.89 -2.29 0.81 1.00 1474 1285
## z1[27] -0.49 -0.46 0.73 0.71 -1.75 0.65 1.00 1168 1054
## z1[28] -0.69 -0.67 0.91 0.87 -2.21 0.81 1.00 1488 1286
## z1[29] -1.69 -1.63 1.04 1.00 -3.48 -0.13 1.00 998 891
## z1[30] -0.39 -0.37 0.77 0.77 -1.70 0.84 1.00 1217 1143
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1
## 0 13 4
## 1 6 7cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1
## 0 9 0
## 1 4 0
##
## , , = b
##
##
## 0 1
## 0 4 4
## 1 2 7par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1
## 0 13 4
## 1 6 7cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1
## 0 9 0
## 1 4 0
##
## , , = b
##
## map
## 0 1
## 0 4 4
## 1 2 7par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
multi logistic regression
# of samples is N,
# of trials of a sample i is mi,
objective variable [0,n], integer
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~B(mi, pi), # of acutual incident
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-3.stan
multi logistic regression
data{
int N;
int K;
array[N] int m;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~binomial_logit(m,z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
}
}n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)
mdl=cmdstan_model('./ex6-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -148.389
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -95.2055 0.00130982 0.00101506 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -95.21
## b[1] -0.23
## b[2] 0.83
## b[3] 1.20
## y1[1] 1.00
## y1[2] 2.00
## y1[3] 2.00
## y1[4] 3.00
## y1[5] 1.00
## y1[6] 7.00
## y1[7] 2.00
## y1[8] 5.00
## y1[9] 1.00
## y1[10] 1.00
## y1[11] 4.00
## y1[12] 2.00
## y1[13] 3.00
## y1[14] 6.00
## y1[15] 2.00
## y1[16] 3.00
## y1[17] 1.00
## y1[18] 4.00
## y1[19] 3.00
## y1[20] 1.00
## y1[21] 0.00
## y1[22] 5.00
## y1[23] 5.00
## y1[24] 3.00
## y1[25] 1.00
## y1[26] 5.00
## y1[27] 3.00
## y1[28] 6.00
## y1[29] 1.00
## y1[30] 5.00
## z1[1] 0.43
## z1[2] 1.57
## z1[3] -0.95
## z1[4] 1.01
## z1[5] 0.20
## z1[6] 0.99
## z1[7] 0.35
## z1[8] -0.14
## z1[9] 1.03
## z1[10] 1.66
## z1[11] 0.04
## z1[12] 0.35
## z1[13] -0.92
## z1[14] 1.19
## z1[15] 0.34
## z1[16] -0.23
## z1[17] 1.57
## z1[18] -0.38
## z1[19] -1.00
## z1[20] -0.44
## z1[21] -0.86
## z1[22] 0.51
## z1[23] 1.49
## z1[24] -0.99
## z1[25] -0.63
## z1[26] 1.22
## z1[27] -0.07
## z1[28] 0.41
## z1[29] 0.53
## z1[30] 0.58mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -96.69 -96.38 1.23 0.97 -99.18 -95.39 1.00 891 1092
## b[1] -0.25 -0.24 0.25 0.26 -0.68 0.16 1.00 1447 1449
## b[2] 0.84 0.84 0.35 0.35 0.29 1.42 1.00 1147 1135
## b[3] 1.24 1.25 0.33 0.34 0.70 1.79 1.00 1148 1229
## y1[1] 1.80 2.00 0.87 1.48 0.00 3.00 1.00 1932 NA
## y1[2] 1.65 2.00 0.55 0.00 1.00 2.00 1.00 2012 NA
## y1[3] 2.51 2.00 1.47 1.48 0.00 5.00 1.00 1843 1976
## y1[4] 2.95 3.00 0.89 1.48 1.00 4.00 1.00 1989 NA
## y1[5] 1.11 1.00 0.70 1.48 0.00 2.00 1.00 2093 NA
## y1[6] 6.58 7.00 1.42 1.48 4.00 9.00 1.00 2008 NA
## y1[7] 3.53 4.00 1.29 1.48 1.00 6.00 1.00 1948 NA
## y1[8] 2.74 3.00 1.26 1.48 1.00 5.00 1.00 1986 2026
## y1[9] 3.66 4.00 1.02 1.48 2.00 5.00 1.00 1824 NA
## y1[10] 0.84 1.00 0.37 0.00 0.00 1.00 1.00 2113 NA
## y1[11] 4.03 4.00 1.55 1.48 2.00 7.00 1.00 1786 1956
## y1[12] 1.16 1.00 0.72 1.48 0.00 2.00 1.00 1924 NA
## y1[13] 1.38 1.00 1.04 1.48 0.00 3.00 1.00 2178 2014
## y1[14] 6.87 7.00 1.34 1.48 5.00 9.00 1.00 1854 NA
## y1[15] 2.89 3.00 1.18 1.48 1.00 5.00 1.00 1754 NA
## y1[16] 3.54 4.00 1.46 1.48 1.00 6.00 1.00 1912 2029
## y1[17] 0.83 1.00 0.38 0.00 0.00 1.00 1.00 2110 NA
## y1[18] 1.97 2.00 1.15 1.48 0.00 4.00 1.00 1882 1900
## y1[19] 1.87 2.00 1.23 1.48 0.00 4.00 1.00 1740 1763
## y1[20] 1.59 2.00 0.98 1.48 0.00 3.00 1.00 2039 1671
## y1[21] 0.30 0.00 0.46 0.00 0.00 1.00 1.00 1736 NA
## y1[22] 4.38 4.00 1.34 1.48 2.00 6.00 1.00 1919 2058
## y1[23] 7.33 8.00 1.23 1.48 5.00 9.00 1.00 1810 NA
## y1[24] 2.20 2.00 1.33 1.48 0.00 4.00 1.00 1717 1801
## y1[25] 0.66 1.00 0.68 1.48 0.00 2.00 1.00 2209 NA
## y1[26] 3.89 4.00 0.97 1.48 2.00 5.00 1.00 1897 NA
## y1[27] 1.93 2.00 1.02 1.48 0.00 4.00 1.00 1944 NA
## y1[28] 4.88 5.00 1.50 1.48 2.00 7.00 1.00 1861 1940
## y1[29] 0.65 1.00 0.48 0.00 0.00 1.00 1.00 1992 NA
## y1[30] 5.13 5.00 1.40 1.48 3.00 7.00 1.00 2012 1939
## z1[1] 0.44 0.45 0.29 0.28 -0.04 0.91 1.00 1702 1534
## z1[2] 1.60 1.58 0.40 0.40 0.95 2.26 1.00 1230 1327
## z1[3] -0.98 -0.96 0.33 0.33 -1.53 -0.45 1.00 1231 1410
## z1[4] 1.03 1.03 0.25 0.25 0.61 1.45 1.00 1658 1653
## z1[5] 0.22 0.22 0.35 0.34 -0.35 0.79 1.00 1565 1569
## z1[6] 1.02 1.02 0.25 0.25 0.60 1.44 1.00 1674 1692
## z1[7] 0.37 0.37 0.31 0.31 -0.13 0.87 1.00 1689 1622
## z1[8] -0.16 -0.16 0.27 0.27 -0.61 0.27 1.00 1423 1448
## z1[9] 1.05 1.05 0.26 0.25 0.63 1.48 1.00 1630 1772
## z1[10] 1.69 1.68 0.43 0.43 1.00 2.41 1.00 1206 1307
## z1[11] 0.02 0.03 0.31 0.31 -0.50 0.54 1.00 1371 1354
## z1[12] 0.34 0.34 0.40 0.40 -0.32 1.00 1.00 1290 1358
## z1[13] -0.95 -0.93 0.32 0.32 -1.49 -0.43 1.00 1238 1401
## z1[14] 1.21 1.21 0.29 0.28 0.75 1.70 1.00 1442 1470
## z1[15] 0.33 0.33 0.40 0.40 -0.33 0.98 1.00 1292 1358
## z1[16] -0.25 -0.24 0.25 0.26 -0.67 0.17 1.00 1448 1449
## z1[17] 1.60 1.59 0.40 0.40 0.95 2.26 1.00 1229 1327
## z1[18] -0.40 -0.40 0.24 0.24 -0.80 -0.01 1.00 1448 1546
## z1[19] -1.03 -1.01 0.35 0.35 -1.60 -0.48 1.00 1219 1378
## z1[20] -0.46 -0.46 0.24 0.24 -0.86 -0.07 1.00 1436 1475
## z1[21] -0.88 -0.87 0.31 0.31 -1.40 -0.40 1.00 1259 1444
## z1[22] 0.53 0.52 0.27 0.27 0.08 0.96 1.00 1768 1626
## z1[23] 1.52 1.51 0.38 0.37 0.92 2.15 1.00 1255 1424
## z1[24] -1.02 -1.01 0.34 0.35 -1.59 -0.47 1.00 1220 1378
## z1[25] -0.66 -0.65 0.26 0.26 -1.10 -0.25 1.00 1356 1467
## z1[26] 1.25 1.24 0.30 0.29 0.77 1.74 1.00 1413 1447
## z1[27] -0.08 -0.08 0.28 0.29 -0.57 0.37 1.00 1402 1536
## z1[28] 0.43 0.43 0.29 0.29 -0.05 0.91 1.00 1692 1585
## z1[29] 0.55 0.55 0.27 0.26 0.11 0.98 1.00 1787 1742
## z1[30] 0.60 0.60 0.26 0.26 0.17 1.01 1.00 1822 1759
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4 5 7 8
## 0 1 2 0 0 0 0 0 0
## 1 0 4 3 0 0 0 0 0
## 2 0 0 4 0 0 0 0 0
## 3 0 1 1 3 1 1 0 0
## 4 0 0 0 0 1 0 0 0
## 5 0 0 0 0 4 1 1 0
## 6 0 0 0 0 0 0 1 0
## 9 0 0 0 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4 5 7 8
## 0 1 1 0 0 0 0 0 0
## 1 0 1 3 0 0 0 0 0
## 2 0 0 3 0 0 0 0 0
## 3 0 1 0 2 0 0 0 0
## 4 0 0 0 0 0 0 0 0
## 5 0 0 0 0 2 0 0 0
## 6 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4 5 7 8
## 0 0 1 0 0 0 0 0 0
## 1 0 3 0 0 0 0 0 0
## 2 0 0 1 0 0 0 0 0
## 3 0 0 1 1 1 1 0 0
## 4 0 0 0 0 1 0 0 0
## 5 0 0 0 0 2 1 1 0
## 6 0 0 0 0 0 0 1 0
## 9 0 0 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 4 5 7 8
## 0 1 2 0 0 0 0 0 0
## 1 1 4 2 0 0 0 0 0
## 2 0 0 4 0 0 0 0 0
## 3 0 1 1 3 1 1 0 0
## 4 0 0 0 0 0 1 0 0
## 5 0 0 0 0 4 1 1 0
## 6 0 0 0 0 0 0 1 0
## 9 0 0 0 0 0 0 0 1cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 4 5 7 8
## 0 1 1 0 0 0 0 0 0
## 1 1 1 2 0 0 0 0 0
## 2 0 0 3 0 0 0 0 0
## 3 0 1 0 2 0 0 0 0
## 4 0 0 0 0 0 0 0 0
## 5 0 0 0 0 2 0 0 0
## 6 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 4 5 7 8
## 0 0 1 0 0 0 0 0 0
## 1 0 3 0 0 0 0 0 0
## 2 0 0 1 0 0 0 0 0
## 3 0 0 1 1 1 1 0 0
## 4 0 0 0 0 0 1 0 0
## 5 0 0 0 0 2 1 1 0
## 6 0 0 0 0 0 0 1 0
## 9 0 0 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
interaction of variable
n=50
mdl=cmdstan_model('./ex5.stan')
tb=tibble(x=runif(n,-3,3),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,x,y,col=ca),
qplot(data=tb,x,y,col=cb),ncol=2)
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)
fn(f0)## Initial log joint probability = -966.17
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -25.9661 0.000395456 0.00220611 1 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -27.95 -27.62 1.45 1.21 -30.84 -26.32 1.00 812 1070
## b[1] 0.35 0.36 0.21 0.21 0.01 0.67 1.00 598 538
## b[2] 0.88 0.87 0.09 0.09 0.72 1.02 1.00 2570 1543
## b[3] 0.33 0.32 0.29 0.28 -0.12 0.82 1.00 566 601
## s 1.08 1.07 0.12 0.11 0.91 1.29 1.00 2142 1403
## y1[1] 0.70 0.70 1.12 1.09 -1.15 2.50 1.00 2036 1943
## y1[2] 1.08 1.08 1.11 1.08 -0.77 2.92 1.00 1968 2026
## y1[3] -2.03 -2.01 1.11 1.12 -3.84 -0.22 1.00 1620 1751
## y1[4] 1.25 1.24 1.12 1.08 -0.58 3.06 1.00 1934 2050
## y1[5] 2.57 2.57 1.13 1.11 0.72 4.50 1.00 1818 1881
## y1[6] 1.97 1.96 1.09 1.07 0.17 3.83 1.00 2038 1900
## y1[7] 0.28 0.30 1.11 1.10 -1.50 2.12 1.00 1838 1855
## y1[8] 2.32 2.31 1.12 1.14 0.53 4.19 1.00 1966 1933
## y1[9] -1.40 -1.40 1.14 1.12 -3.28 0.46 1.00 1814 1914
## y1[10] -2.28 -2.30 1.16 1.15 -4.18 -0.35 1.00 1818 1959
## y1[11] 1.58 1.60 1.14 1.15 -0.32 3.44 1.00 1896 1882
## y1[12] 0.13 0.14 1.11 1.07 -1.69 1.96 1.00 1990 1952
## y1[13] -0.89 -0.89 1.09 1.05 -2.67 0.93 1.00 1870 1899
## y1[14] 0.23 0.22 1.11 1.09 -1.56 2.08 1.00 2098 1933
## y1[15] 2.24 2.26 1.14 1.16 0.41 4.06 1.00 1913 2103
## y1[16] -0.03 -0.03 1.11 1.11 -1.83 1.88 1.00 1843 1822
## y1[17] -0.83 -0.81 1.09 1.05 -2.64 0.92 1.00 2094 2032
## y1[18] 0.23 0.23 1.11 1.09 -1.54 2.07 1.00 1981 1892
## y1[19] 0.56 0.56 1.11 1.11 -1.24 2.46 1.00 1918 1773
## y1[20] 1.69 1.69 1.10 1.09 -0.15 3.55 1.00 1740 1884
## y1[21] -0.29 -0.28 1.08 1.06 -2.08 1.48 1.00 1947 1857
## y1[22] 2.97 2.98 1.15 1.11 1.10 4.87 1.00 1897 1822
## y1[23] -0.58 -0.57 1.11 1.12 -2.42 1.23 1.00 1916 2016
## y1[24] 2.99 3.01 1.16 1.12 1.08 4.93 1.00 2089 1818
## y1[25] -1.38 -1.36 1.13 1.09 -3.29 0.42 1.00 1940 1958
## y1[26] 1.84 1.85 1.12 1.10 0.00 3.67 1.00 1965 1703
## y1[27] 1.62 1.62 1.12 1.09 -0.26 3.43 1.00 1843 1893
## y1[28] -1.64 -1.62 1.12 1.08 -3.49 0.19 1.00 1826 1916
## y1[29] 2.93 2.93 1.13 1.08 1.10 4.77 1.00 1985 1786
## y1[30] -1.65 -1.62 1.10 1.10 -3.54 0.09 1.00 1772 1840
## y1[31] 1.31 1.28 1.10 1.07 -0.49 3.11 1.00 1869 1861
## y1[32] 0.35 0.37 1.10 1.09 -1.50 2.18 1.00 1834 1961
## y1[33] -0.75 -0.76 1.11 1.08 -2.60 1.06 1.00 2027 1910
## y1[34] -0.74 -0.75 1.09 1.05 -2.52 1.07 1.00 1976 1972
## y1[35] -1.53 -1.53 1.10 1.10 -3.29 0.28 1.00 1925 1863
## y1[36] -0.08 -0.10 1.09 1.07 -1.92 1.73 1.00 1936 2002
## y1[37] -0.42 -0.40 1.10 1.09 -2.26 1.33 1.00 1888 1879
## y1[38] 1.11 1.14 1.09 1.07 -0.64 2.92 1.00 1920 1891
## y1[39] -0.55 -0.56 1.11 1.08 -2.35 1.28 1.00 2117 1848
## y1[40] -1.12 -1.11 1.10 1.05 -2.97 0.69 1.00 1688 1934
## y1[41] -2.25 -2.27 1.12 1.13 -4.06 -0.44 1.00 1784 1821
## y1[42] 0.31 0.31 1.12 1.07 -1.54 2.18 1.00 1921 1658
## y1[43] 1.62 1.60 1.09 1.04 -0.15 3.44 1.00 1891 1866
## y1[44] 2.48 2.47 1.14 1.16 0.65 4.40 1.00 1845 1521
## y1[45] 1.36 1.39 1.11 1.12 -0.49 3.08 1.00 2037 1893
## y1[46] -1.48 -1.48 1.14 1.11 -3.32 0.37 1.00 2132 1813
## y1[47] -1.13 -1.14 1.11 1.08 -2.93 0.72 1.00 1973 1969
## y1[48] -1.96 -1.98 1.14 1.13 -3.80 -0.14 1.00 1782 1893
## y1[49] -0.61 -0.60 1.09 1.09 -2.43 1.17 1.00 1826 1965
## y1[50] 0.78 0.81 1.10 1.05 -1.09 2.57 1.00 1735 1776
## m1[1] 0.73 0.73 0.21 0.21 0.39 1.07 1.00 1259 1333
## m1[2] 1.05 1.05 0.22 0.21 0.69 1.40 1.00 1300 1286
## m1[3] -1.98 -1.98 0.29 0.28 -2.46 -1.49 1.00 1358 1357
## m1[4] 1.26 1.27 0.24 0.24 0.86 1.65 1.00 631 759
## m1[5] 2.60 2.59 0.34 0.33 2.03 3.15 1.00 849 1120
## m1[6] 2.01 2.01 0.26 0.26 1.59 2.45 1.00 1535 1224
## m1[7] 0.28 0.28 0.21 0.21 -0.08 0.63 1.00 1262 1279
## m1[8] 2.33 2.33 0.28 0.28 1.87 2.80 1.00 1620 1220
## m1[9] -1.36 -1.36 0.29 0.29 -1.84 -0.90 1.00 1752 1445
## m1[10] -2.24 -2.24 0.31 0.30 -2.75 -1.73 1.00 1515 1581
## m1[11] 1.56 1.57 0.26 0.26 1.13 1.98 1.00 665 978
## m1[12] 0.13 0.14 0.21 0.21 -0.20 0.45 1.00 611 546
## m1[13] -0.92 -0.91 0.26 0.26 -1.35 -0.49 1.00 1511 1446
## m1[14] 0.23 0.23 0.21 0.21 -0.13 0.57 1.00 1263 1279
## m1[15] 2.28 2.28 0.32 0.30 1.76 2.80 1.00 787 990
## m1[16] -0.01 -0.01 0.22 0.21 -0.37 0.35 1.00 1286 1436
## m1[17] -0.86 -0.87 0.22 0.22 -1.23 -0.49 1.00 822 955
## m1[18] 0.19 0.20 0.21 0.21 -0.14 0.51 1.00 607 502
## m1[19] 0.55 0.56 0.21 0.21 0.20 0.88 1.00 595 627
## m1[20] 1.69 1.69 0.27 0.26 1.24 2.13 1.00 685 957
## m1[21] -0.29 -0.29 0.23 0.22 -0.67 0.09 1.00 1305 1368
## m1[22] 2.99 2.99 0.33 0.32 2.44 3.55 1.00 1787 1184
## m1[23] -0.59 -0.59 0.21 0.21 -0.94 -0.24 1.00 737 728
## m1[24] 2.99 2.99 0.33 0.32 2.44 3.54 1.00 1786 1184
## m1[25] -1.40 -1.41 0.25 0.25 -1.82 -0.98 1.00 1052 1282
## m1[26] 1.82 1.83 0.25 0.25 1.41 2.25 1.00 1492 1247
## m1[27] 1.60 1.61 0.26 0.26 1.16 2.03 1.00 671 894
## m1[28] -1.64 -1.64 0.31 0.32 -2.14 -1.13 1.00 1835 1435
## m1[29] 2.92 2.92 0.33 0.32 2.38 3.47 1.00 1774 1198
## m1[30] -1.69 -1.70 0.27 0.26 -2.14 -1.24 1.00 1204 1282
## m1[31] 1.31 1.31 0.25 0.24 0.90 1.70 1.00 635 815
## m1[32] 0.34 0.34 0.21 0.21 -0.02 0.68 1.00 1255 1222
## m1[33] -0.74 -0.74 0.25 0.25 -1.15 -0.33 1.00 1431 1403
## m1[34] -0.72 -0.72 0.22 0.22 -1.07 -0.36 1.00 776 844
## m1[35] -1.49 -1.48 0.30 0.30 -1.97 -1.00 1.00 1791 1413
## m1[36] -0.08 -0.07 0.20 0.21 -0.41 0.24 1.00 635 567
## m1[37] -0.41 -0.41 0.23 0.23 -0.80 -0.03 1.00 1331 1369
## m1[38] 1.07 1.07 0.22 0.21 0.71 1.42 1.00 1303 1286
## m1[39] -0.57 -0.57 0.24 0.24 -0.96 -0.18 1.00 1369 1330
## m1[40] -1.13 -1.14 0.24 0.23 -1.53 -0.74 1.00 929 1075
## m1[41] -2.25 -2.25 0.31 0.30 -2.76 -1.74 1.00 1520 1556
## m1[42] 0.32 0.32 0.21 0.21 -0.04 0.66 1.00 1256 1222
## m1[43] 1.60 1.61 0.24 0.23 1.21 2.00 1.00 1423 1196
## m1[44] 2.47 2.47 0.33 0.31 1.92 3.01 1.00 824 1081
## m1[45] 1.40 1.40 0.23 0.22 1.02 1.77 1.00 1368 1347
## m1[46] -1.44 -1.45 0.25 0.25 -1.87 -1.01 1.00 1069 1282
## m1[47] -1.11 -1.10 0.27 0.27 -1.56 -0.67 1.00 1617 1499
## m1[48] -1.95 -1.95 0.29 0.27 -2.43 -1.47 1.00 1341 1320
## m1[49] -0.62 -0.62 0.21 0.21 -0.97 -0.27 1.00 747 728
## m1[50] 0.81 0.82 0.22 0.22 0.45 1.15 1.00 601 665
fn(f1)## Initial log joint probability = -60.4972
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -25.2382 0.000130537 0.0018028 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -27.78 -27.47 1.62 1.51 -30.75 -25.82 1.00 840 1101
## b[1] 0.13 0.14 0.27 0.28 -0.31 0.56 1.01 708 1084
## b[2] 0.88 0.88 0.09 0.09 0.74 1.03 1.00 2780 1882
## b[3] 0.46 0.45 0.32 0.31 -0.06 0.97 1.01 693 1001
## b[4] 0.37 0.38 0.32 0.32 -0.16 0.90 1.01 820 916
## s 1.07 1.07 0.12 0.11 0.90 1.28 1.00 2431 1653
## y1[1] 0.63 0.63 1.09 1.05 -1.22 2.40 1.00 2009 1931
## y1[2] 1.38 1.38 1.14 1.11 -0.49 3.24 1.00 1937 1971
## y1[3] -1.83 -1.81 1.12 1.08 -3.75 -0.02 1.00 1979 1729
## y1[4] 1.44 1.47 1.16 1.12 -0.50 3.30 1.00 1923 1750
## y1[5] 2.79 2.81 1.11 1.11 1.01 4.59 1.00 1696 2023
## y1[6] 1.92 1.90 1.12 1.06 0.06 3.80 1.00 2029 1854
## y1[7] 0.21 0.20 1.11 1.10 -1.64 2.04 1.00 1961 1969
## y1[8] 2.29 2.28 1.15 1.15 0.42 4.18 1.00 1786 1748
## y1[9] -1.45 -1.44 1.14 1.13 -3.36 0.45 1.00 1875 1880
## y1[10] -2.10 -2.11 1.12 1.08 -4.01 -0.34 1.00 2036 1877
## y1[11] 1.37 1.37 1.13 1.07 -0.53 3.25 1.00 1698 1927
## y1[12] 0.29 0.27 1.08 1.07 -1.45 2.10 1.00 2120 2011
## y1[13] -1.03 -1.04 1.10 1.09 -2.86 0.75 1.00 1730 1745
## y1[14] 0.50 0.49 1.12 1.08 -1.32 2.39 1.00 1901 1955
## y1[15] 2.07 2.04 1.17 1.15 0.15 3.92 1.00 1567 1807
## y1[16] -0.10 -0.09 1.10 1.09 -1.92 1.71 1.00 1819 1966
## y1[17] -1.10 -1.10 1.12 1.13 -2.92 0.74 1.00 1860 1858
## y1[18] 0.37 0.32 1.11 1.11 -1.44 2.23 1.00 1527 1776
## y1[19] 0.71 0.70 1.09 1.10 -1.07 2.49 1.00 1908 1878
## y1[20] 1.46 1.46 1.11 1.08 -0.38 3.28 1.00 1842 1767
## y1[21] -0.05 -0.07 1.14 1.08 -1.95 1.84 1.00 1854 1693
## y1[22] 2.96 2.95 1.15 1.11 1.12 4.79 1.00 1719 1971
## y1[23] -0.43 -0.45 1.12 1.09 -2.25 1.43 1.00 1829 1933
## y1[24] 2.90 2.90 1.13 1.12 1.01 4.77 1.00 1927 1952
## y1[25] -1.25 -1.21 1.14 1.08 -3.16 0.59 1.00 1882 1790
## y1[26] 2.12 2.13 1.16 1.13 0.20 4.00 1.00 1921 1929
## y1[27] 1.40 1.41 1.11 1.15 -0.44 3.19 1.00 1742 1881
## y1[28] -1.77 -1.74 1.11 1.11 -3.57 0.06 1.00 1672 1863
## y1[29] 2.87 2.85 1.16 1.17 0.95 4.76 1.00 2035 1787
## y1[30] -1.92 -1.95 1.13 1.13 -3.75 -0.05 1.00 1723 1955
## y1[31] 1.48 1.48 1.12 1.07 -0.35 3.33 1.00 1830 1831
## y1[32] 0.21 0.22 1.12 1.13 -1.59 2.09 1.00 2025 1936
## y1[33] -0.87 -0.86 1.10 1.09 -2.65 0.94 1.00 2032 1976
## y1[34] -0.62 -0.61 1.08 1.07 -2.41 1.14 1.00 2168 2004
## y1[35] -1.61 -1.62 1.14 1.15 -3.44 0.41 1.00 2018 2035
## y1[36] -0.31 -0.30 1.10 1.08 -2.14 1.49 1.00 2014 2015
## y1[37] -0.14 -0.16 1.15 1.10 -2.06 1.78 1.00 1937 1821
## y1[38] 0.96 0.97 1.13 1.10 -0.92 2.82 1.00 1963 1895
## y1[39] -0.65 -0.67 1.11 1.06 -2.45 1.24 1.00 1757 1881
## y1[40] -1.42 -1.41 1.11 1.09 -3.26 0.47 1.00 1742 1712
## y1[41] -2.10 -2.10 1.12 1.10 -3.93 -0.25 1.00 1932 1899
## y1[42] 0.58 0.56 1.12 1.08 -1.24 2.47 1.00 1995 1713
## y1[43] 1.50 1.48 1.11 1.14 -0.34 3.34 1.00 1978 2013
## y1[44] 2.64 2.64 1.10 1.11 0.81 4.43 1.00 1975 1822
## y1[45] 1.31 1.33 1.10 1.09 -0.47 3.06 1.00 1936 1893
## y1[46] -1.30 -1.29 1.09 1.09 -3.10 0.41 1.01 1582 1574
## y1[47] -1.23 -1.25 1.12 1.13 -3.04 0.62 1.00 2106 1903
## y1[48] -2.17 -2.15 1.17 1.13 -4.07 -0.31 1.00 1921 1972
## y1[49] -0.46 -0.45 1.12 1.10 -2.36 1.41 1.00 2031 1867
## y1[50] 0.62 0.64 1.12 1.10 -1.20 2.47 1.00 1918 2019
## m1[1] 0.64 0.64 0.23 0.23 0.26 1.04 1.00 1681 1603
## m1[2] 1.34 1.35 0.33 0.32 0.78 1.87 1.01 894 1005
## m1[3] -1.84 -1.84 0.32 0.32 -2.37 -1.32 1.00 1896 1802
## m1[4] 1.43 1.42 0.29 0.30 0.94 1.90 1.00 1484 1569
## m1[5] 2.77 2.76 0.38 0.39 2.16 3.39 1.00 1754 1712
## m1[6] 1.93 1.93 0.28 0.27 1.49 2.39 1.00 2168 1535
## m1[7] 0.19 0.19 0.23 0.23 -0.19 0.58 1.00 1596 1606
## m1[8] 2.26 2.26 0.29 0.29 1.78 2.74 1.00 2295 1753
## m1[9] -1.47 -1.47 0.31 0.30 -1.98 -0.98 1.00 1694 1749
## m1[10] -2.11 -2.11 0.34 0.33 -2.67 -1.56 1.00 1988 1852
## m1[11] 1.35 1.35 0.31 0.30 0.85 1.84 1.00 868 1310
## m1[12] 0.29 0.29 0.26 0.25 -0.14 0.71 1.00 1398 1294
## m1[13] -1.02 -1.02 0.28 0.28 -1.48 -0.55 1.00 1629 1780
## m1[14] 0.51 0.52 0.32 0.31 -0.03 1.03 1.01 873 1062
## m1[15] 2.08 2.08 0.35 0.35 1.51 2.65 1.00 1048 1542
## m1[16] -0.11 -0.10 0.24 0.24 -0.50 0.29 1.00 1572 1537
## m1[17] -1.09 -1.09 0.29 0.29 -1.55 -0.61 1.01 760 1062
## m1[18] 0.35 0.35 0.26 0.25 -0.09 0.77 1.00 1398 1293
## m1[19] 0.71 0.71 0.27 0.26 0.26 1.14 1.00 1406 1299
## m1[20] 1.48 1.48 0.31 0.31 0.97 1.98 1.00 897 1367
## m1[21] -0.01 0.00 0.32 0.32 -0.56 0.51 1.01 912 1132
## m1[22] 2.92 2.93 0.34 0.35 2.37 3.48 1.00 2503 1819
## m1[23] -0.44 -0.43 0.26 0.26 -0.86 -0.02 1.00 1483 1530
## m1[24] 2.92 2.92 0.34 0.35 2.36 3.47 1.00 2502 1819
## m1[25] -1.26 -1.25 0.29 0.28 -1.73 -0.79 1.00 1700 1750
## m1[26] 2.12 2.12 0.36 0.34 1.53 2.70 1.01 994 1146
## m1[27] 1.39 1.39 0.31 0.31 0.88 1.88 1.00 877 1329
## m1[28] -1.75 -1.74 0.33 0.32 -2.29 -1.23 1.00 1741 1814
## m1[29] 2.85 2.86 0.34 0.34 2.31 3.40 1.00 2484 1820
## m1[30] -1.92 -1.92 0.33 0.33 -2.46 -1.38 1.01 902 1126
## m1[31] 1.47 1.47 0.30 0.30 0.98 1.95 1.00 1493 1568
## m1[32] 0.24 0.24 0.23 0.23 -0.14 0.63 1.00 1604 1606
## m1[33] -0.84 -0.84 0.27 0.26 -1.29 -0.39 1.00 1615 1703
## m1[34] -0.57 -0.57 0.27 0.26 -1.00 -0.14 1.00 1509 1601
## m1[35] -1.59 -1.59 0.31 0.31 -2.11 -1.09 1.00 1713 1816
## m1[36] -0.30 -0.29 0.27 0.27 -0.74 0.14 1.01 702 1073
## m1[37] -0.14 -0.12 0.33 0.32 -0.69 0.39 1.01 925 1217
## m1[38] 0.98 0.98 0.24 0.23 0.60 1.39 1.00 1796 1603
## m1[39] -0.67 -0.67 0.26 0.26 -1.11 -0.24 1.00 1593 1637
## m1[40] -1.36 -1.36 0.30 0.30 -1.84 -0.85 1.01 803 1118
## m1[41] -2.11 -2.11 0.34 0.33 -2.67 -1.57 1.00 1991 1852
## m1[42] 0.60 0.61 0.32 0.31 0.06 1.13 1.01 869 1072
## m1[43] 1.52 1.53 0.26 0.24 1.10 1.97 1.00 1992 1651
## m1[44] 2.64 2.63 0.37 0.38 2.04 3.24 1.00 1728 1713
## m1[45] 1.31 1.32 0.25 0.24 0.91 1.74 1.00 1932 1536
## m1[46] -1.30 -1.29 0.29 0.28 -1.78 -0.83 1.00 1712 1719
## m1[47] -1.21 -1.21 0.29 0.28 -1.70 -0.73 1.00 1661 1749
## m1[48] -2.18 -2.18 0.35 0.35 -2.75 -1.61 1.01 960 1211
## m1[49] -0.47 -0.47 0.26 0.26 -0.90 -0.05 1.00 1488 1530
## m1[50] 0.60 0.60 0.28 0.28 0.16 1.03 1.00 743 1166
fn(f2)## Initial log joint probability = -285.114
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 37 -23.3564 6.67346e-05 0.000830854 0.3772 0.9909 46
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.61 -26.29 1.81 1.72 -30.07 -24.24 1.00 687 1057
## b[1] -0.18 -0.18 0.33 0.31 -0.71 0.36 1.00 635 996
## b[2] 0.89 0.89 0.09 0.09 0.75 1.03 1.00 2164 1312
## b[3] 0.94 0.94 0.41 0.41 0.27 1.61 1.00 529 796
## b[4] 0.89 0.90 0.43 0.40 0.20 1.62 1.00 561 935
## b[5] -1.21 -1.21 0.64 0.66 -2.21 -0.16 1.01 414 563
## s 1.05 1.04 0.12 0.12 0.88 1.26 1.00 2322 1499
## y1[1] 0.85 0.85 1.09 1.07 -0.96 2.65 1.00 1936 1813
## y1[2] 0.81 0.82 1.16 1.12 -1.08 2.74 1.00 1647 1819
## y1[3] -1.65 -1.68 1.08 1.10 -3.46 0.08 1.00 1970 2042
## y1[4] 1.64 1.63 1.11 1.08 -0.18 3.48 1.00 2018 1833
## y1[5] 3.02 3.02 1.09 1.06 1.26 4.84 1.00 1900 1860
## y1[6] 2.13 2.14 1.10 1.10 0.32 3.90 1.00 2106 1542
## y1[7] 0.41 0.41 1.07 1.04 -1.28 2.16 1.00 1960 2014
## y1[8] 2.43 2.45 1.08 1.05 0.64 4.20 1.00 1983 1930
## y1[9] -1.33 -1.36 1.07 1.06 -3.12 0.47 1.00 1869 1932
## y1[10] -1.95 -1.95 1.10 1.10 -3.70 -0.11 1.00 1791 1453
## y1[11] 1.02 1.03 1.10 1.11 -0.75 2.88 1.00 1822 1733
## y1[12] 0.49 0.49 1.09 1.10 -1.26 2.26 1.00 1731 1794
## y1[13] -0.89 -0.89 1.11 1.10 -2.72 0.92 1.00 1759 1621
## y1[14] -0.01 -0.01 1.13 1.11 -1.88 1.84 1.00 1663 1903
## y1[15] 1.81 1.83 1.14 1.15 -0.11 3.69 1.00 2007 1847
## y1[16] 0.04 0.04 1.06 1.07 -1.73 1.75 1.00 1945 1879
## y1[17] -1.46 -1.43 1.10 1.11 -3.33 0.30 1.00 1769 1793
## y1[18] 0.55 0.52 1.09 1.11 -1.21 2.35 1.00 2039 1954
## y1[19] 0.92 0.93 1.10 1.12 -0.84 2.67 1.00 1883 1973
## y1[20] 1.21 1.19 1.12 1.12 -0.58 3.03 1.00 1627 1777
## y1[21] -0.56 -0.48 1.13 1.08 -2.41 1.24 1.00 1948 1750
## y1[22] 3.14 3.11 1.11 1.08 1.36 4.93 1.00 1963 2010
## y1[23] -0.21 -0.23 1.07 1.06 -1.97 1.59 1.00 1750 1643
## y1[24] 3.14 3.16 1.10 1.10 1.35 4.99 1.00 1915 1746
## y1[25] -1.05 -1.04 1.10 1.06 -2.92 0.74 1.00 1870 1892
## y1[26] 1.66 1.68 1.14 1.10 -0.26 3.53 1.00 1791 2012
## y1[27] 1.09 1.08 1.10 1.11 -0.72 2.88 1.00 1545 1708
## y1[28] -1.63 -1.65 1.10 1.08 -3.37 0.18 1.00 2063 1845
## y1[29] 3.05 3.08 1.13 1.11 1.20 4.83 1.00 1945 1984
## y1[30] -2.27 -2.23 1.11 1.11 -4.11 -0.46 1.00 1642 1668
## y1[31] 1.66 1.67 1.12 1.09 -0.14 3.48 1.00 1855 1738
## y1[32] 0.42 0.40 1.07 1.10 -1.29 2.15 1.00 2092 2055
## y1[33] -0.70 -0.69 1.13 1.14 -2.53 1.11 1.00 2057 1918
## y1[34] -0.38 -0.39 1.09 1.05 -2.13 1.41 1.00 2115 1904
## y1[35] -1.44 -1.41 1.12 1.11 -3.28 0.37 1.00 1966 1867
## y1[36] -0.60 -0.63 1.11 1.09 -2.38 1.26 1.00 1606 1919
## y1[37] -0.66 -0.67 1.15 1.12 -2.48 1.27 1.00 2062 1899
## y1[38] 1.19 1.19 1.09 1.06 -0.58 2.95 1.00 2183 2027
## y1[39] -0.52 -0.52 1.09 1.07 -2.29 1.28 1.00 2024 1855
## y1[40] -1.69 -1.68 1.13 1.14 -3.59 0.12 1.00 1528 1913
## y1[41] -1.91 -1.92 1.10 1.05 -3.71 -0.06 1.00 1979 1972
## y1[42] 0.07 0.08 1.15 1.11 -1.85 2.01 1.00 1498 1631
## y1[43] 1.77 1.77 1.07 1.07 0.01 3.53 1.00 1967 1974
## y1[44] 2.92 2.91 1.12 1.08 1.12 4.70 1.00 1828 2013
## y1[45] 1.46 1.44 1.07 1.05 -0.30 3.24 1.00 1860 1792
## y1[46] -1.12 -1.10 1.09 1.05 -2.90 0.70 1.00 1943 1739
## y1[47] -1.11 -1.11 1.10 1.08 -2.93 0.69 1.00 1996 1889
## y1[48] -2.53 -2.53 1.12 1.07 -4.41 -0.67 1.00 1584 1666
## y1[49] -0.31 -0.29 1.08 1.09 -2.10 1.45 1.00 1634 1719
## y1[50] 0.30 0.31 1.10 1.11 -1.54 2.08 1.00 1834 2109
## m1[1] 0.82 0.82 0.25 0.24 0.40 1.23 1.00 1488 1411
## m1[2] 0.83 0.83 0.44 0.42 0.14 1.55 1.00 1115 1072
## m1[3] -1.65 -1.66 0.33 0.34 -2.16 -1.10 1.00 1645 1161
## m1[4] 1.64 1.65 0.32 0.33 1.11 2.16 1.00 1355 1497
## m1[5] 3.01 3.01 0.40 0.40 2.33 3.65 1.00 1514 1627
## m1[6] 2.12 2.12 0.29 0.28 1.65 2.60 1.00 1572 1550
## m1[7] 0.36 0.36 0.25 0.24 -0.05 0.77 1.00 1494 1402
## m1[8] 2.45 2.45 0.31 0.30 1.94 2.95 1.00 1615 1452
## m1[9] -1.32 -1.32 0.31 0.31 -1.81 -0.79 1.00 1737 1441
## m1[10] -1.92 -1.93 0.34 0.35 -2.46 -1.35 1.00 1713 1280
## m1[11] 1.05 1.06 0.35 0.34 0.47 1.64 1.00 754 1116
## m1[12] 0.50 0.50 0.28 0.29 0.03 0.96 1.00 1308 1245
## m1[13] -0.86 -0.86 0.29 0.29 -1.31 -0.37 1.00 1647 1510
## m1[14] -0.01 -0.01 0.43 0.42 -0.71 0.73 1.00 1052 1158
## m1[15] 1.79 1.80 0.39 0.37 1.17 2.46 1.00 885 1278
## m1[16] 0.06 0.06 0.25 0.25 -0.36 0.48 1.00 1513 1423
## m1[17] -1.41 -1.42 0.35 0.33 -1.96 -0.85 1.00 651 1133
## m1[18] 0.56 0.56 0.28 0.29 0.08 1.01 1.00 1309 1246
## m1[19] 0.92 0.92 0.29 0.29 0.44 1.39 1.00 1312 1214
## m1[20] 1.19 1.20 0.36 0.35 0.60 1.79 1.00 776 1116
## m1[21] -0.54 -0.54 0.44 0.42 -1.26 0.22 1.00 1028 1219
## m1[22] 3.12 3.12 0.35 0.34 2.56 3.70 1.00 1707 1373
## m1[23] -0.24 -0.24 0.28 0.28 -0.69 0.23 1.00 1356 1302
## m1[24] 3.12 3.12 0.35 0.34 2.56 3.69 1.00 1707 1373
## m1[25] -1.07 -1.07 0.30 0.31 -1.55 -0.56 1.00 1504 1253
## m1[26] 1.62 1.61 0.45 0.45 0.89 2.36 1.00 1207 1194
## m1[27] 1.09 1.10 0.36 0.34 0.51 1.68 1.00 760 1117
## m1[28] -1.60 -1.61 0.33 0.33 -2.12 -1.04 1.00 1791 1394
## m1[29] 3.05 3.05 0.35 0.33 2.49 3.62 1.00 1698 1389
## m1[30] -2.25 -2.26 0.38 0.36 -2.85 -1.61 1.00 726 1170
## m1[31] 1.69 1.70 0.32 0.33 1.15 2.21 1.00 1359 1490
## m1[32] 0.41 0.41 0.25 0.24 0.00 0.82 1.00 1494 1362
## m1[33] -0.68 -0.68 0.28 0.28 -1.13 -0.21 1.00 1613 1487
## m1[34] -0.37 -0.37 0.28 0.28 -0.82 0.10 1.00 1377 1159
## m1[35] -1.44 -1.45 0.32 0.31 -1.95 -0.90 1.00 1766 1394
## m1[36] -0.61 -0.61 0.33 0.31 -1.14 -0.07 1.00 627 1001
## m1[37] -0.66 -0.67 0.44 0.42 -1.39 0.10 1.00 1026 1235
## m1[38] 1.16 1.16 0.26 0.25 0.74 1.57 1.00 1489 1456
## m1[39] -0.51 -0.51 0.27 0.27 -0.94 -0.06 1.00 1581 1406
## m1[40] -1.69 -1.70 0.35 0.34 -2.25 -1.10 1.00 670 1135
## m1[41] -1.93 -1.94 0.34 0.35 -2.47 -1.36 1.00 1715 1280
## m1[42] 0.08 0.08 0.43 0.42 -0.62 0.82 1.00 1057 1143
## m1[43] 1.71 1.71 0.27 0.26 1.25 2.15 1.00 1513 1444
## m1[44] 2.87 2.88 0.39 0.39 2.21 3.50 1.00 1497 1656
## m1[45] 1.50 1.50 0.26 0.25 1.06 1.94 1.00 1499 1367
## m1[46] -1.11 -1.11 0.30 0.31 -1.59 -0.59 1.00 1514 1271
## m1[47] -1.06 -1.06 0.30 0.30 -1.52 -0.55 1.00 1687 1511
## m1[48] -2.52 -2.52 0.39 0.38 -3.14 -1.85 1.00 756 1150
## m1[49] -0.27 -0.28 0.28 0.29 -0.72 0.20 1.00 1361 1219
## m1[50] 0.29 0.30 0.33 0.32 -0.26 0.83 1.00 665 988
tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,xa,y,col=ca),
qplot(data=tb,xa,y,col=cb),
qplot(data=tb,xb,y,col=ca),
qplot(data=tb,xb,y,col=cb))
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)
fn(f0)## Initial log joint probability = -83.506
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -49.0092 0.000177649 0.00121943 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -51.86 -51.46 1.96 1.71 -55.60 -49.40 1.00 774 1028
## b[1] 1.03 1.03 0.41 0.39 0.35 1.69 1.01 866 902
## b[2] 1.03 1.03 0.26 0.25 0.61 1.46 1.00 2134 1495
## b[3] 1.01 1.00 0.25 0.24 0.60 1.41 1.00 2136 1270
## b[4] 0.04 0.05 0.50 0.46 -0.78 0.86 1.00 929 939
## b[5] 2.07 2.06 0.52 0.51 1.23 2.94 1.01 832 879
## s 1.76 1.74 0.20 0.19 1.47 2.11 1.01 1913 1489
## y1[1] 0.74 0.75 1.89 1.84 -2.36 3.87 1.00 1876 1971
## y1[2] -0.33 -0.29 1.84 1.86 -3.33 2.60 1.00 1699 1841
## y1[3] 0.51 0.49 1.86 1.80 -2.54 3.48 1.00 1818 1954
## y1[4] 4.98 5.07 1.84 1.82 1.92 7.90 1.00 1857 1986
## y1[5] 1.56 1.54 1.88 1.86 -1.60 4.80 1.00 1878 1998
## y1[6] 3.09 3.06 1.83 1.79 0.05 6.11 1.00 1892 1931
## y1[7] 1.25 1.29 1.83 1.81 -1.79 4.23 1.00 1874 1797
## y1[8] 1.24 1.26 1.88 1.88 -1.85 4.30 1.00 1978 2058
## y1[9] 1.91 1.90 1.81 1.79 -0.96 4.92 1.00 1927 1701
## y1[10] 4.64 4.68 1.92 1.91 1.52 7.73 1.00 1880 1808
## y1[11] 1.04 1.05 1.80 1.76 -1.96 3.97 1.00 1866 1858
## y1[12] 0.61 0.65 1.86 1.77 -2.54 3.65 1.00 2091 2030
## y1[13] 2.70 2.71 1.82 1.79 -0.24 5.61 1.00 1562 1968
## y1[14] 2.56 2.59 1.86 1.82 -0.44 5.60 1.00 1907 1873
## y1[15] 1.51 1.50 1.82 1.74 -1.47 4.48 1.00 2001 1933
## y1[16] 4.12 4.18 1.86 1.85 1.05 7.16 1.00 1924 1805
## y1[17] -1.23 -1.23 1.87 1.82 -4.33 1.80 1.00 1737 1768
## y1[18] 0.49 0.51 1.80 1.73 -2.59 3.30 1.00 1906 1796
## y1[19] 6.00 6.01 2.01 1.97 2.63 9.29 1.00 1919 1933
## y1[20] 2.70 2.76 1.85 1.84 -0.45 5.70 1.00 1945 1776
## y1[21] 1.80 1.85 1.85 1.80 -1.37 4.75 1.00 2026 1877
## y1[22] 2.45 2.41 1.90 1.88 -0.60 5.55 1.00 1909 1846
## y1[23] 2.17 2.15 1.87 1.89 -0.84 5.17 1.00 1909 1755
## y1[24] -0.45 -0.49 1.83 1.79 -3.43 2.54 1.00 1869 1972
## y1[25] 0.59 0.59 1.90 1.89 -2.50 3.63 1.00 1852 1857
## y1[26] 2.64 2.68 1.84 1.75 -0.56 5.64 1.00 1901 1846
## y1[27] 3.83 3.84 1.90 1.89 0.73 6.92 1.00 1791 1900
## y1[28] 2.79 2.84 1.86 1.90 -0.30 5.76 1.00 1758 1756
## y1[29] -0.21 -0.22 1.83 1.79 -3.27 2.84 1.00 1832 1973
## y1[30] 0.65 0.65 1.85 1.80 -2.41 3.71 1.00 2010 2101
## y1[31] 3.14 3.16 1.82 1.78 0.24 6.11 1.00 2091 1880
## y1[32] 2.03 2.00 1.83 1.82 -1.04 4.99 1.00 1756 1956
## y1[33] 1.49 1.56 1.81 1.72 -1.57 4.36 1.00 1650 1777
## y1[34] 1.66 1.64 1.80 1.76 -1.34 4.60 1.00 1809 1969
## y1[35] 3.54 3.56 1.91 1.88 0.40 6.70 1.00 1952 1890
## y1[36] 2.95 2.97 1.90 1.88 -0.18 6.08 1.00 1915 2056
## y1[37] 2.05 2.02 1.84 1.78 -1.05 5.01 1.00 2040 1900
## y1[38] -0.87 -0.88 1.80 1.75 -3.92 2.08 1.00 2019 1929
## y1[39] 4.32 4.29 1.87 1.83 1.24 7.49 1.00 1982 2058
## y1[40] 1.97 1.94 1.78 1.74 -0.89 4.89 1.00 1974 1956
## y1[41] -2.65 -2.62 1.93 1.91 -5.83 0.55 1.00 1719 2015
## y1[42] 0.66 0.64 1.86 1.85 -2.39 3.75 1.00 1691 1596
## y1[43] 2.90 2.94 1.84 1.83 -0.10 5.93 1.00 1921 1971
## y1[44] -1.89 -1.91 1.90 1.87 -4.96 1.25 1.00 2076 1917
## y1[45] 2.74 2.72 1.83 1.72 -0.32 5.75 1.00 1738 1888
## y1[46] 1.14 1.10 1.86 1.85 -1.80 4.24 1.00 1862 1928
## y1[47] 1.41 1.43 1.84 1.79 -1.64 4.37 1.00 1885 1793
## y1[48] -0.48 -0.45 1.84 1.78 -3.47 2.44 1.00 2040 1732
## y1[49] 3.69 3.60 1.97 1.97 0.54 6.94 1.00 1644 1800
## y1[50] 1.60 1.55 1.82 1.84 -1.35 4.56 1.00 1862 1713
## m1[1] 0.78 0.78 0.54 0.55 -0.10 1.64 1.00 1196 1199
## m1[2] -0.34 -0.33 0.45 0.44 -1.08 0.41 1.01 948 903
## m1[3] 0.50 0.49 0.53 0.52 -0.38 1.38 1.00 1289 1187
## m1[4] 5.02 5.02 0.63 0.62 4.00 6.06 1.00 1346 1086
## m1[5] 1.56 1.56 0.53 0.51 0.71 2.42 1.00 1321 1287
## m1[6] 3.05 3.05 0.51 0.50 2.23 3.87 1.00 1445 1356
## m1[7] 1.25 1.25 0.57 0.57 0.33 2.18 1.00 1466 1404
## m1[8] 1.31 1.31 0.65 0.63 0.18 2.36 1.00 1331 1375
## m1[9] 1.90 1.90 0.53 0.53 1.03 2.76 1.00 1345 1453
## m1[10] 4.63 4.65 0.64 0.63 3.60 5.67 1.00 1747 1476
## m1[11] 1.04 1.02 0.49 0.49 0.24 1.84 1.00 1194 1317
## m1[12] 0.56 0.56 0.58 0.58 -0.35 1.48 1.00 1689 1592
## m1[13] 2.75 2.73 0.46 0.45 2.01 3.52 1.00 1184 1384
## m1[14] 2.56 2.54 0.53 0.52 1.70 3.40 1.00 1303 1400
## m1[15] 1.51 1.51 0.55 0.55 0.61 2.41 1.00 1272 1459
## m1[16] 4.09 4.09 0.56 0.54 3.19 5.02 1.00 1338 1230
## m1[17] -1.21 -1.22 0.59 0.57 -2.21 -0.25 1.00 1228 1063
## m1[18] 0.49 0.49 0.59 0.59 -0.44 1.43 1.00 1745 1568
## m1[19] 5.98 5.96 0.80 0.79 4.68 7.28 1.00 1511 1284
## m1[20] 2.68 2.69 0.59 0.58 1.73 3.60 1.00 1298 1052
## m1[21] 1.73 1.73 0.61 0.62 0.71 2.71 1.00 1471 1451
## m1[22] 2.43 2.42 0.61 0.63 1.45 3.42 1.01 1574 1450
## m1[23] 2.25 2.26 0.51 0.50 1.43 3.07 1.00 1163 1131
## m1[24] -0.48 -0.49 0.49 0.48 -1.28 0.38 1.01 1121 1009
## m1[25] 0.56 0.55 0.66 0.67 -0.49 1.64 1.00 1423 1441
## m1[26] 2.62 2.61 0.50 0.51 1.81 3.45 1.00 1384 1485
## m1[27] 3.84 3.83 0.55 0.54 2.94 4.75 1.00 1372 1299
## m1[28] 2.87 2.86 0.46 0.45 2.13 3.62 1.00 1318 1454
## m1[29] -0.24 -0.25 0.50 0.49 -1.02 0.61 1.00 1348 1508
## m1[30] 0.66 0.66 0.49 0.48 -0.14 1.47 1.00 1189 1224
## m1[31] 3.14 3.13 0.55 0.56 2.25 4.02 1.01 2060 1424
## m1[32] 1.99 2.00 0.51 0.49 1.15 2.80 1.00 1089 1193
## m1[33] 1.53 1.53 0.43 0.43 0.82 2.23 1.01 966 1039
## m1[34] 1.67 1.68 0.46 0.45 0.91 2.43 1.00 1013 1251
## m1[35] 3.53 3.54 0.70 0.71 2.42 4.67 1.01 2423 1486
## m1[36] 2.96 2.96 0.56 0.57 2.05 3.86 1.00 2024 1485
## m1[37] 2.05 2.04 0.62 0.62 1.02 3.06 1.00 1884 1524
## m1[38] -0.85 -0.86 0.55 0.54 -1.74 0.09 1.00 1534 1485
## m1[39] 4.43 4.45 0.52 0.52 3.57 5.27 1.00 1775 1560
## m1[40] 1.95 1.96 0.51 0.49 1.10 2.78 1.00 1103 1194
## m1[41] -2.55 -2.57 0.78 0.77 -3.83 -1.27 1.00 1812 1554
## m1[42] 0.68 0.68 0.47 0.48 -0.07 1.45 1.00 1143 1324
## m1[43] 2.92 2.92 0.56 0.56 2.03 3.83 1.00 1513 1422
## m1[44] -1.88 -1.88 0.69 0.67 -3.03 -0.73 1.00 1696 1380
## m1[45] 2.72 2.72 0.44 0.43 2.01 3.43 1.00 1427 1387
## m1[46] 1.15 1.13 0.63 0.63 0.14 2.16 1.00 1500 1135
## m1[47] 1.43 1.43 0.54 0.55 0.55 2.32 1.00 1419 1273
## m1[48] -0.43 -0.43 0.50 0.49 -1.28 0.36 1.00 1036 993
## m1[49] 3.65 3.65 0.78 0.79 2.38 4.87 1.00 1253 1366
## m1[50] 1.53 1.51 0.63 0.63 0.52 2.56 1.00 1542 1370
fn(f1)## Initial log joint probability = -477.205
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 28 -46.7072 0.000155759 0.00087093 1 1 31
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -50.15 -49.80 2.10 1.98 -54.11 -47.46 1.00 688 964
## b[1] 0.67 0.67 0.45 0.43 -0.06 1.40 1.00 714 810
## b[2] 1.11 1.12 0.25 0.26 0.70 1.53 1.00 1891 1691
## b[3] 1.08 1.08 0.25 0.23 0.66 1.48 1.00 1934 1280
## b[4] 1.00 0.99 0.71 0.71 -0.19 2.17 1.01 505 652
## b[5] 2.89 2.89 0.64 0.65 1.86 3.89 1.00 570 778
## b[6] -1.99 -1.94 1.03 1.03 -3.73 -0.38 1.00 411 527
## s 1.70 1.68 0.19 0.18 1.42 2.04 1.00 1621 1556
## y1[1] 1.33 1.31 1.80 1.82 -1.60 4.23 1.00 1716 1865
## y1[2] -0.84 -0.86 1.79 1.74 -3.75 2.24 1.00 1580 1930
## y1[3] 0.09 0.08 1.82 1.87 -2.85 2.99 1.00 2004 1775
## y1[4] 4.63 4.63 1.82 1.77 1.72 7.60 1.00 1900 1900
## y1[5] 1.89 1.90 1.75 1.75 -0.98 4.82 1.00 1876 1786
## y1[6] 3.49 3.51 1.75 1.73 0.63 6.27 1.00 1853 1969
## y1[7] 0.55 0.57 1.81 1.85 -2.44 3.53 1.00 1831 1970
## y1[8] 0.98 1.00 1.80 1.78 -1.96 3.94 1.00 1985 1857
## y1[9] 2.28 2.26 1.81 1.75 -0.61 5.22 1.00 2018 1891
## y1[10] 5.29 5.31 1.86 1.83 2.28 8.32 1.00 1805 1932
## y1[11] 0.67 0.70 1.79 1.79 -2.31 3.58 1.00 1828 1685
## y1[12] 1.09 1.11 1.79 1.80 -1.90 3.91 1.00 1729 1931
## y1[13] 2.13 2.09 1.81 1.74 -0.87 5.13 1.00 2030 2015
## y1[14] 1.93 1.96 1.80 1.72 -1.07 4.77 1.00 1896 1894
## y1[15] 2.18 2.19 1.76 1.72 -0.65 5.05 1.00 1906 1604
## y1[16] 3.54 3.52 1.79 1.68 0.59 6.54 1.00 2008 1785
## y1[17] -1.80 -1.76 1.78 1.81 -4.74 1.15 1.00 1772 2037
## y1[18] 0.94 0.90 1.79 1.76 -2.00 3.93 1.00 1807 1921
## y1[19] 5.65 5.68 1.96 1.91 2.28 8.82 1.00 1867 1957
## y1[20] 2.45 2.43 1.78 1.79 -0.34 5.44 1.00 1819 1822
## y1[21] 2.06 2.02 1.82 1.77 -0.87 5.04 1.00 1807 1973
## y1[22] 2.79 2.86 1.80 1.81 -0.22 5.67 1.00 2054 1697
## y1[23] 2.01 1.99 1.76 1.77 -0.86 4.88 1.00 1970 1762
## y1[24] -0.91 -0.88 1.77 1.75 -3.77 1.95 1.00 1625 1855
## y1[25] 1.11 1.14 1.85 1.81 -2.05 4.13 1.00 1728 1917
## y1[26] 3.06 3.06 1.76 1.79 0.14 5.87 1.00 1841 2036
## y1[27] 3.26 3.30 1.83 1.83 0.30 6.20 1.00 1979 1856
## y1[28] 3.37 3.34 1.78 1.65 0.45 6.31 1.00 1809 1856
## y1[29] 0.14 0.16 1.79 1.76 -2.89 3.09 1.00 1805 1889
## y1[30] 0.29 0.29 1.74 1.69 -2.58 3.12 1.00 1962 1891
## y1[31] 3.52 3.54 1.78 1.70 0.52 6.43 1.00 1940 1930
## y1[32] 1.71 1.67 1.76 1.75 -1.25 4.58 1.00 1909 1973
## y1[33] 1.19 1.22 1.75 1.74 -1.70 3.99 1.00 1985 1981
## y1[34] 1.31 1.30 1.71 1.63 -1.45 4.14 1.00 1720 2105
## y1[35] 4.07 4.08 1.85 1.87 1.08 7.06 1.00 2037 1973
## y1[36] 3.46 3.42 1.86 1.86 0.43 6.55 1.00 1502 1420
## y1[37] 2.39 2.37 1.87 1.80 -0.67 5.47 1.00 1802 1930
## y1[38] -0.39 -0.39 1.80 1.80 -3.32 2.46 1.00 1706 1706
## y1[39] 4.98 4.97 1.78 1.69 2.06 7.98 1.00 1908 1883
## y1[40] 1.64 1.60 1.78 1.74 -1.30 4.56 1.00 1794 1924
## y1[41] -2.12 -2.08 1.91 1.93 -5.25 1.10 1.00 1813 1823
## y1[42] 1.24 1.23 1.80 1.70 -1.78 4.25 1.00 1878 1721
## y1[43] 2.27 2.32 1.79 1.78 -0.65 5.27 1.00 1850 1839
## y1[44] -1.55 -1.55 1.76 1.74 -4.41 1.40 1.00 2098 1843
## y1[45] 3.13 3.14 1.79 1.77 0.14 6.10 1.00 1972 1893
## y1[46] 0.38 0.42 1.84 1.90 -2.69 3.39 1.00 1728 1840
## y1[47] 0.79 0.78 1.82 1.83 -2.20 3.88 1.00 2096 1904
## y1[48] -0.92 -0.97 1.80 1.78 -3.80 2.06 1.00 1711 1827
## y1[49] 4.41 4.46 1.94 1.91 1.29 7.68 1.00 1624 1882
## y1[50] 1.21 1.22 1.80 1.80 -1.72 4.12 1.00 1924 1815
## m1[1] 1.36 1.35 0.60 0.59 0.35 2.37 1.00 898 1187
## m1[2] -0.81 -0.82 0.51 0.49 -1.64 0.05 1.00 661 947
## m1[3] 0.10 0.09 0.56 0.56 -0.80 1.01 1.00 934 1314
## m1[4] 4.59 4.60 0.66 0.64 3.50 5.65 1.00 1925 1486
## m1[5] 1.91 1.91 0.54 0.52 1.02 2.76 1.00 1668 1592
## m1[6] 3.50 3.50 0.52 0.51 2.63 4.33 1.00 1065 1337
## m1[7] 0.53 0.54 0.66 0.64 -0.55 1.62 1.00 1080 1138
## m1[8] 0.96 0.94 0.66 0.66 -0.10 2.03 1.00 1200 1230
## m1[9] 2.26 2.27 0.54 0.54 1.38 3.14 1.00 1294 1543
## m1[10] 5.20 5.21 0.67 0.66 4.11 6.32 1.00 1192 1463
## m1[11] 0.68 0.68 0.52 0.51 -0.16 1.53 1.00 904 1191
## m1[12] 1.11 1.10 0.60 0.59 0.12 2.12 1.00 1192 1445
## m1[13] 2.15 2.15 0.55 0.53 1.25 3.05 1.01 1131 1439
## m1[14] 1.94 1.93 0.60 0.59 0.95 2.91 1.00 1196 1444
## m1[15] 2.13 2.14 0.61 0.61 1.12 3.17 1.00 931 1177
## m1[16] 3.59 3.59 0.60 0.58 2.58 4.54 1.00 1536 1354
## m1[17] -1.75 -1.76 0.65 0.61 -2.83 -0.66 1.00 771 1133
## m1[18] 1.04 1.03 0.61 0.61 0.03 2.07 1.00 1242 1462
## m1[19] 5.62 5.63 0.80 0.77 4.32 6.93 1.01 2349 1624
## m1[20] 2.45 2.45 0.58 0.57 1.49 3.45 1.00 1263 1218
## m1[21] 2.07 2.07 0.60 0.62 1.08 3.05 1.00 1430 1573
## m1[22] 2.83 2.83 0.61 0.60 1.81 3.81 1.00 1335 1673
## m1[23] 1.98 1.98 0.52 0.50 1.14 2.86 1.00 1085 1137
## m1[24] -0.95 -0.96 0.54 0.53 -1.82 -0.02 1.00 730 1097
## m1[25] 1.12 1.11 0.70 0.69 -0.03 2.29 1.00 1099 1380
## m1[26] 3.04 3.04 0.52 0.51 2.19 3.87 1.00 1104 1388
## m1[27] 3.31 3.33 0.60 0.58 2.30 4.28 1.00 1485 1222
## m1[28] 3.30 3.30 0.48 0.46 2.50 4.06 1.00 1002 1182
## m1[29] 0.25 0.25 0.53 0.51 -0.63 1.11 1.00 1139 1287
## m1[30] 0.28 0.28 0.52 0.52 -0.57 1.14 1.00 863 1228
## m1[31] 3.60 3.58 0.58 0.59 2.67 4.57 1.00 1472 1313
## m1[32] 1.69 1.69 0.52 0.52 0.87 2.58 1.01 1046 1298
## m1[33] 1.20 1.20 0.46 0.46 0.45 1.95 1.00 821 1036
## m1[34] 1.36 1.35 0.49 0.49 0.59 2.17 1.01 901 1085
## m1[35] 4.03 4.02 0.72 0.71 2.86 5.22 1.00 1679 1512
## m1[36] 3.41 3.39 0.58 0.59 2.47 4.38 1.00 1536 1352
## m1[37] 2.43 2.42 0.63 0.64 1.38 3.46 1.00 1816 1562
## m1[38] -0.40 -0.41 0.56 0.55 -1.33 0.50 1.00 1398 1323
## m1[39] 4.99 4.98 0.57 0.57 4.04 5.93 1.00 1036 1222
## m1[40] 1.66 1.66 0.53 0.52 0.82 2.56 1.01 1051 1335
## m1[41] -2.24 -2.25 0.75 0.76 -3.48 -1.03 1.00 2032 1261
## m1[42] 1.24 1.23 0.53 0.51 0.39 2.12 1.00 865 1047
## m1[43] 2.32 2.34 0.62 0.60 1.26 3.31 1.00 1318 1133
## m1[44] -1.51 -1.51 0.69 0.69 -2.66 -0.40 1.00 1711 1303
## m1[45] 3.15 3.14 0.47 0.47 2.38 3.94 1.00 1130 1402
## m1[46] 0.42 0.43 0.71 0.70 -0.71 1.58 1.00 1185 1434
## m1[47] 0.72 0.73 0.64 0.62 -0.32 1.79 1.00 1076 1130
## m1[48] -0.91 -0.91 0.55 0.53 -1.82 0.01 1.00 715 1128
## m1[49] 4.44 4.44 0.86 0.87 3.04 5.86 1.00 902 1073
## m1[50] 1.21 1.21 0.63 0.63 0.18 2.24 1.00 1209 1555
fn(f2)## Initial log joint probability = -123.462
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 35 -22.216 6.34509e-05 0.00153826 0.2246 1 42
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.78 -26.42 2.24 2.12 -30.98 -23.80 1.00 773 1160
## b[1] 0.23 0.23 0.28 0.27 -0.23 0.70 1.00 980 1464
## b[2] 0.96 0.95 0.15 0.15 0.71 1.21 1.00 2509 1773
## b[3] 0.84 0.84 0.15 0.14 0.59 1.09 1.00 2380 1292
## b[4] 1.46 1.45 0.44 0.44 0.75 2.18 1.00 756 953
## b[5] 2.20 2.20 0.41 0.41 1.55 2.87 1.00 618 914
## b[6] -1.15 -1.15 0.14 0.14 -1.38 -0.91 1.00 2011 1770
## b[7] -1.13 -1.14 0.65 0.68 -2.23 -0.04 1.00 483 612
## s 1.05 1.04 0.12 0.12 0.87 1.27 1.00 2022 1496
## y1[1] 1.90 1.91 1.14 1.12 0.03 3.77 1.00 1962 2014
## y1[2] -1.49 -1.46 1.12 1.08 -3.35 0.31 1.00 1924 1931
## y1[3] 1.15 1.14 1.09 1.10 -0.64 2.97 1.00 1902 2012
## y1[4] 3.45 3.46 1.16 1.18 1.52 5.33 1.00 1816 1894
## y1[5] 0.65 0.66 1.15 1.16 -1.27 2.50 1.00 2010 2058
## y1[6] 3.35 3.36 1.11 1.11 1.53 5.16 1.00 1967 1997
## y1[7] 0.25 0.24 1.12 1.10 -1.57 2.04 1.00 1628 1792
## y1[8] 3.54 3.56 1.19 1.20 1.63 5.45 1.00 1974 1846
## y1[9] 1.73 1.71 1.11 1.12 -0.15 3.53 1.00 1964 1881
## y1[10] 4.69 4.67 1.14 1.09 2.82 6.59 1.00 1889 1911
## y1[11] 1.12 1.10 1.10 1.11 -0.76 2.97 1.00 1563 1845
## y1[12] 3.06 3.07 1.16 1.09 1.07 5.06 1.00 2131 2012
## y1[13] 2.39 2.40 1.09 1.08 0.64 4.17 1.00 1927 1974
## y1[14] 2.60 2.56 1.16 1.10 0.73 4.52 1.00 1974 1971
## y1[15] 2.95 2.94 1.14 1.09 1.07 4.78 1.00 1569 1495
## y1[16] 4.04 4.03 1.12 1.05 2.15 5.91 1.00 1810 1800
## y1[17] -2.54 -2.55 1.13 1.10 -4.39 -0.69 1.00 1735 1894
## y1[18] 3.14 3.15 1.18 1.16 1.19 5.06 1.00 2034 1885
## y1[19] 3.37 3.34 1.18 1.16 1.43 5.38 1.00 2286 2056
## y1[20] 1.27 1.27 1.11 1.07 -0.52 3.06 1.00 2125 1872
## y1[21] 2.32 2.31 1.17 1.15 0.38 4.21 1.00 1827 1871
## y1[22] 3.95 3.98 1.15 1.14 2.01 5.79 1.00 2097 1660
## y1[23] 0.90 0.88 1.08 1.06 -0.88 2.66 1.00 1920 2046
## y1[24] -1.41 -1.42 1.13 1.15 -3.25 0.44 1.00 1778 1955
## y1[25] 3.17 3.20 1.15 1.13 1.26 5.07 1.00 1733 1626
## y1[26] 2.80 2.77 1.10 1.06 1.07 4.62 1.00 1885 1700
## y1[27] 4.18 4.16 1.12 1.11 2.36 6.04 1.00 1979 2030
## y1[28] 2.54 2.56 1.09 1.06 0.73 4.31 1.00 1868 1477
## y1[29] 0.09 0.11 1.07 1.01 -1.69 1.87 1.00 1978 2054
## y1[30] 0.76 0.74 1.10 1.04 -1.03 2.61 1.00 1971 1937
## y1[31] 4.05 4.01 1.13 1.13 2.16 5.94 1.00 1967 1973
## y1[32] 1.23 1.22 1.12 1.11 -0.56 3.07 1.00 1832 1757
## y1[33] 0.63 0.61 1.11 1.09 -1.21 2.51 1.00 1978 1971
## y1[34] 0.89 0.89 1.10 1.06 -0.91 2.71 1.00 1966 1876
## y1[35] 6.46 6.44 1.18 1.13 4.51 8.36 1.00 2022 1908
## y1[36] 3.98 3.97 1.09 1.04 2.19 5.74 1.00 2018 1970
## y1[37] 3.36 3.36 1.12 1.05 1.44 5.20 1.00 1953 2011
## y1[38] -0.91 -0.93 1.11 1.09 -2.72 0.89 1.00 2048 2035
## y1[39] 3.35 3.37 1.15 1.10 1.40 5.22 1.00 1715 1911
## y1[40] 1.30 1.29 1.11 1.10 -0.46 3.09 1.00 1720 1910
## y1[41] -5.17 -5.19 1.22 1.18 -7.12 -3.16 1.00 1987 1967
## y1[42] 1.36 1.37 1.11 1.11 -0.43 3.14 1.00 1824 2014
## y1[43] 4.13 4.11 1.17 1.20 2.22 6.08 1.00 1784 1835
## y1[44] -3.24 -3.24 1.16 1.17 -5.13 -1.34 1.00 2029 1981
## y1[45] 2.25 2.28 1.12 1.07 0.37 4.06 1.00 1755 1704
## y1[46] 0.40 0.41 1.16 1.11 -1.48 2.30 1.00 2091 1849
## y1[47] 0.53 0.51 1.15 1.14 -1.35 2.37 1.00 1679 1787
## y1[48] -1.15 -1.13 1.09 1.09 -2.97 0.57 1.00 1862 2051
## y1[49] 2.11 2.12 1.18 1.18 0.22 4.02 1.00 1996 1827
## y1[50] 3.12 3.11 1.15 1.12 1.19 5.00 1.00 1659 1887
## m1[1] 1.89 1.88 0.36 0.38 1.30 2.48 1.00 1334 1392
## m1[2] -1.51 -1.50 0.33 0.32 -2.05 -0.95 1.00 1048 1497
## m1[3] 1.12 1.11 0.37 0.36 0.53 1.74 1.00 1073 1495
## m1[4] 3.44 3.44 0.43 0.43 2.75 4.15 1.00 2613 1607
## m1[5] 0.63 0.64 0.36 0.36 0.05 1.22 1.00 1272 1439
## m1[6] 3.34 3.34 0.34 0.33 2.78 3.90 1.00 1366 1567
## m1[7] 0.24 0.24 0.41 0.41 -0.42 0.93 1.00 1337 1595
## m1[8] 3.53 3.52 0.53 0.50 2.66 4.40 1.00 1435 1341
## m1[9] 1.69 1.69 0.35 0.34 1.11 2.26 1.00 1392 1460
## m1[10] 4.72 4.73 0.43 0.42 4.02 5.41 1.00 1397 1531
## m1[11] 1.13 1.13 0.32 0.31 0.62 1.68 1.00 1035 1411
## m1[12] 3.05 3.05 0.45 0.42 2.32 3.79 1.00 1800 1579
## m1[13] 2.37 2.37 0.35 0.35 1.79 2.93 1.00 1355 1407
## m1[14] 2.63 2.62 0.39 0.40 2.01 3.25 1.00 1380 1422
## m1[15] 2.97 2.96 0.39 0.39 2.34 3.60 1.00 1376 1394
## m1[16] 4.03 4.03 0.39 0.38 3.39 4.67 1.00 1916 1528
## m1[17] -2.56 -2.56 0.42 0.42 -3.27 -1.88 1.00 1240 1350
## m1[18] 3.18 3.17 0.47 0.44 2.42 3.95 1.00 1861 1603
## m1[19] 3.33 3.33 0.57 0.57 2.42 4.26 1.00 2865 1670
## m1[20] 1.26 1.26 0.37 0.36 0.65 1.87 1.00 1808 1207
## m1[21] 2.31 2.31 0.38 0.38 1.68 2.93 1.00 1936 1627
## m1[22] 3.96 3.97 0.41 0.41 3.27 4.61 1.00 2410 1676
## m1[23] 0.91 0.91 0.33 0.33 0.38 1.45 1.00 1458 1238
## m1[24] -1.44 -1.44 0.34 0.34 -2.01 -0.87 1.00 1107 1584
## m1[25] 3.14 3.13 0.49 0.49 2.36 3.96 1.00 1693 1617
## m1[26] 2.80 2.81 0.33 0.33 2.26 3.34 1.00 1370 1522
## m1[27] 4.20 4.20 0.40 0.40 3.53 4.86 1.00 1727 1619
## m1[28] 2.58 2.57 0.32 0.31 2.04 3.12 1.00 1125 1349
## m1[29] 0.11 0.11 0.33 0.32 -0.43 0.66 1.00 1431 1488
## m1[30] 0.81 0.80 0.33 0.32 0.29 1.36 1.00 1026 1370
## m1[31] 4.07 4.06 0.36 0.35 3.48 4.67 1.00 1688 1585
## m1[32] 1.23 1.23 0.32 0.33 0.71 1.75 1.00 1239 1375
## m1[33] 0.61 0.61 0.29 0.28 0.14 1.08 1.00 1095 1297
## m1[34] 0.90 0.90 0.30 0.30 0.41 1.39 1.00 1122 1354
## m1[35] 6.47 6.48 0.54 0.53 5.60 7.37 1.00 2443 1579
## m1[36] 3.99 3.98 0.36 0.36 3.39 4.60 1.00 1762 1585
## m1[37] 3.35 3.36 0.40 0.39 2.69 3.99 1.00 2167 1536
## m1[38] -0.90 -0.90 0.35 0.35 -1.48 -0.32 1.00 1678 1623
## m1[39] 3.38 3.38 0.41 0.41 2.70 4.03 1.00 1155 1422
## m1[40] 1.31 1.30 0.32 0.32 0.78 1.84 1.00 1220 1388
## m1[41] -5.15 -5.13 0.59 0.59 -6.13 -4.20 1.00 2175 1356
## m1[42] 1.40 1.39 0.33 0.33 0.88 1.95 1.00 1236 1241
## m1[43] 4.13 4.13 0.46 0.46 3.39 4.89 1.00 1469 1255
## m1[44] -3.21 -3.21 0.47 0.48 -4.00 -2.47 1.00 2061 1389
## m1[45] 2.20 2.20 0.31 0.31 1.70 2.71 1.00 1082 1375
## m1[46] 0.37 0.36 0.43 0.43 -0.32 1.10 1.00 1391 1530
## m1[47] 0.54 0.53 0.39 0.39 -0.10 1.20 1.00 1314 1518
## m1[48] -1.14 -1.14 0.35 0.34 -1.71 -0.56 1.00 1052 1222
## m1[49] 2.15 2.14 0.59 0.58 1.17 3.13 1.00 1364 1445
## m1[50] 3.12 3.12 0.46 0.45 2.37 3.90 1.00 1399 1443
